2020.8.11 BRBを再解析(log2FCを計算するため)、それに合わせて再解析
print(Sys.Date())
[1] "2020-08-11"
print(sessionInfo(),locale=FALSE)
R version 4.0.1 (2020-06-06)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux
Matrix products: default
BLAS/LAPACK: /usr/local/intel2018_up1/compilers_and_libraries_2018.0.128/linux/mkl/lib/intel64_lin/libmkl_intel_lp64.so
attached base packages:
[1] grid stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] stringr_1.4.0 hrbrthemes_0.8.0 ggrepel_0.8.2
[4] ggpubr_0.4.0.999 gplots_3.0.4 DESeq2_1.28.1
[7] GGally_2.0.0 vcd_1.4-7 BiocParallel_1.22.0
[10] Matrix_1.2-18 SummarizedExperiment_1.18.2 DelayedArray_0.14.1
[13] matrixStats_0.56.0 motifmatchr_1.10.0 org.Mm.eg.db_3.11.4
[16] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0 org.Hs.eg.db_3.11.4 TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
[19] GenomicFeatures_1.40.1 AnnotationDbi_1.50.3 Biobase_2.48.0
[22] ChIPseeker_1.24.0 clusterProfiler_3.16.0 BSgenome.Mmusculus.UCSC.mm10_1.4.0
[25] ggsignif_0.6.0 chromVAR_1.10.0 purrr_0.3.4
[28] RColorBrewer_1.1-2 ggsci_2.9 readr_1.3.1
[31] tidyr_1.1.1 dplyr_1.0.1 ggplot2_3.3.2
[34] TFBSTools_1.26.0 BSgenome_1.56.0 rtracklayer_1.48.0
[37] Biostrings_2.56.0 XVector_0.28.0 GenomicRanges_1.40.0
[40] GenomeInfoDb_1.24.2 IRanges_2.22.2 S4Vectors_0.26.1
[43] BiocGenerics_0.34.0
loaded via a namespace (and not attached):
[1] rappdirs_0.3.1 R.methodsS3_1.8.0 bit64_4.0.2 knitr_1.29 irlba_2.3.3
[6] R.utils_2.9.2 data.table_1.13.0 KEGGREST_1.28.0 RCurl_1.98-1.2 generics_0.0.2
[11] snow_0.4-3 cowplot_1.0.0 lambda.r_1.2.4 RSQLite_2.2.0 europepmc_0.4
[16] bit_4.0.4 enrichplot_1.8.1 xml2_1.3.2 httpuv_1.5.4 isoband_0.2.2
[21] assertthat_0.2.1 DirichletMultinomial_1.30.0 viridis_0.5.1 xfun_0.16 hms_0.5.3
[26] evaluate_0.14 promises_1.1.1 fansi_0.4.1 progress_1.2.2 caTools_1.18.0
[31] dbplyr_1.4.4 readxl_1.3.1 igraph_1.2.5 DBI_1.1.0 geneplotter_1.66.0
[36] htmlwidgets_1.5.1 futile.logger_1.4.3 reshape_0.8.8 ellipsis_0.3.1 backports_1.1.8
[41] annotate_1.66.0 biomaRt_2.44.1 vctrs_0.3.2 abind_1.4-5 withr_2.2.0
[46] ggforce_0.3.2 triebeard_0.3.0 GenomicAlignments_1.24.0 prettyunits_1.1.1 DOSE_3.14.0
[51] lazyeval_0.2.2 seqLogo_1.54.3 crayon_1.3.4 genefilter_1.70.0 labeling_0.3
[56] pkgconfig_2.0.3 tweenr_1.0.1 nlme_3.1-148 rlang_0.4.7 lifecycle_0.2.0
[61] miniUI_0.1.1.1 downloader_0.4 extrafontdb_1.0 BiocFileCache_1.12.1 cellranger_1.1.0
[66] polyclip_1.10-0 lmtest_0.9-37 urltools_1.7.3 carData_3.0-4 boot_1.3-25
[71] zoo_1.8-8 base64enc_0.1-3 pheatmap_1.0.12 ggridges_0.5.2 png_0.1-7
[76] viridisLite_0.3.0 bitops_1.0-6 R.oo_1.23.0 KernSmooth_2.23-17 blob_1.2.1
[81] qvalue_2.20.0 rstatix_0.6.0 gridGraphics_0.5-0 CNEr_1.24.0 scales_1.1.1
[86] memoise_1.1.0 magrittr_1.5 plyr_1.8.6 gdata_2.18.0 zlibbioc_1.34.0
[91] compiler_4.0.1 scatterpie_0.1.4 plotrix_3.7-8 Rsamtools_2.4.0 cli_2.0.2
[96] formatR_1.7 mgcv_1.8-31 MASS_7.3-51.6 tidyselect_1.1.0 stringi_1.4.6
[101] forcats_0.5.0 yaml_2.2.1 GOSemSim_2.14.1 askpass_1.1 locfit_1.5-9.4
[106] fastmatch_1.1-0 tools_4.0.1 rio_0.5.16 rstudioapi_0.11 TFMPvalue_0.0.8
[111] foreign_0.8-80 gridExtra_2.3 farver_2.0.3 ggraph_2.0.3 digest_0.6.25
[116] rvcheck_0.1.8 BiocManager_1.30.10 FNN_1.1.3 shiny_1.5.0 pracma_2.2.9
[121] Rcpp_1.0.5 car_3.0-8 broom_0.7.0 later_1.1.0.1 gdtools_0.2.2
[126] httr_1.4.2 colorspace_1.4-1 XML_3.99-0.5 splines_4.0.1 uwot_0.1.8
[131] graphlayouts_0.7.0 ggplotify_0.0.5 systemfonts_0.2.3 plotly_4.9.2.1 xtable_1.8-4
[136] jsonlite_1.7.0 futile.options_1.0.1 poweRlaw_0.70.6 tidygraph_1.2.0 R6_2.4.1
[141] pillar_1.4.6 htmltools_0.5.0 mime_0.9 glue_1.4.1 fastmap_1.0.1
[146] DT_0.15 fgsea_1.14.0 utf8_1.1.4 lattice_0.20-41 tibble_3.0.3
[151] curl_4.3 gtools_3.8.2 Rttf2pt1_1.3.8 zip_2.1.0 GO.db_3.11.4
[156] openxlsx_4.1.5 openssl_1.4.2 survival_3.2-3 rmarkdown_2.3 munsell_0.5.0
[161] DO.db_2.9 GenomeInfoDbData_1.2.3 msigdbr_7.1.1 haven_2.3.1 reshape2_1.4.4
[166] gtable_0.3.0 extrafont_0.17
#{r setup 0, include=FALSE} knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)
source("/home/guestA/n70275b/work/rscripts/geomNorm.R")
## ラベルあり
ggpoints <- function(x,...)
ggplot(x,...) + geom_point(stroke=1) +
ggrepel::geom_text_repel(size=4) + theme_minimal() + mycolor
maxchrom <- 19 # 19: mouse, 22: human
mycolor <- ggsci::scale_color_aaas()
# PCA/UMAP
scalerows <- TRUE # gene-wise scaling (pattern is the matter?)
ntop <- 500 # number of top-n genes with high variance
seed <- 123 # set another number if UMAP looks not good
n_nei <- 6 # number of neighboring data points in UMAP #ここをどうしたらいい?
#----------------#
#cluster_num <- 4
filepath_summary <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__TSS_pm5kb_20200624.count.txt"
# deftable 修正版
filepath_def <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/deftable_multicov_ChIL01100111_20200501_3T3_EGFP18_UI_DoxMinus_H3p3K27acK4Kme327me3_withATAC.txt"
#--- サンプルの選択 ---#
folder_name <- "ChIL H3.3, H3K27ac, H3K4me3, H3K27me3, ATAC" #"ChIL_H3K27me3_H3K27me3"
#remove_sample=c("Doxminus_UI_ATAC_4","Doxplus_UI_ATAC_4","Doxminus_D48_ATAC_1","Doxminus_D48_ATAC_4","Doxplus_D48_ATAC_4") #このサンプルを削除(20190917) <ATACの場合>
#use <- quo(!(sample %in% remove_sample)) #このサンプルを削除(20190917)
#--- multiBamSummary の略 ---#
# データ保存用のpath
#csvfilepath <- basename(filepath_summary) %>% sub(".count.txt", "__", .)
csvfilepath <- basename(filepath_summary) %>% sub(".count.txt", "", .) %>% sub("ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__", "", .)
print(csvfilepath)
[1] "TSS_pm5kb_20200624"
filepath_BRBensemble <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/ensemble_list_useast.csv" #BRBの時のgene名リスト
ensemble <- readr::read_csv(filepath_BRBensemble) %>% mutate_if(is.double, as.integer)
Parsed with column specification:
cols(
ens_gene = col_character(),
ext_gene = col_character(),
biotype = col_character(),
chr = col_character()
)
annotate <- partial(right_join,ensemble,by="ens_gene") #2gunで使う
UCSCの形式の場合 (20191016)
# 20200617
def_list <- readr::read_tsv(filepath_def) %>% mutate(seq=factor(seq, c("ATAC","H3p3", "H3K27ac","H3K4me3","H3K27me3"))) %>%
mutate(time1 = time, rep1=rep) %>% unite(time1,rep1,col="time_replicate") %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% mutate(type=factor(type,c("DoxPlus","DoxMinus"))) %>% mutate(rep=factor(rep, c("1", "2", "3", "4")))%>% mutate(time_replicate=factor(time_replicate,c("UI_1", "UI_2", "UI_3", "UI_4", "0h_1","0h_2","24h_1","24h_2","48h_1","48h_2","48h_3","48h_4")))
Parsed with column specification:
cols(
file = col_character(),
multicov_No = col_double(),
sample = col_character(),
group = col_character(),
time = col_character(),
type = col_character(),
seq = col_character(),
rep = col_double()
)
#def_list <- readr::read_tsv(filepath_def) %>% mutate(seq=factor(seq, c("H3p3", "H3K27ac","H3K4me3","H3K27me3"))) %>%
#mutate(time1 = time, rep1=rep) %>% unite(time1,rep1,col="time_replicate") %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% mutate(type=factor(type,c("DoxPlus","DoxMinus"))) %>% mutate(rep=factor(rep, c("1", "2")))%>% mutate(time_replicate=factor(time_replicate,c("UI_1", "UI_2", "0h_1","0h_2","24h_1","24h_2","48h_1","48h_2")))
# add 20200617
groups <- c(
"ATAC_UI_DoxMinus","ATAC_UI_DoxPlus",
"ATAC_48h_DoxMinus","ATAC_48h_DoxPlus",
"H3p3_UI_DoxMinus","H3p3_UI_DoxPlus",
"H3p3_0h_DoxMinus","H3p3_0h_DoxPlus",
"H3p3_24h_DoxMinus","H3p3_24h_DoxPlus",
"H3p3_48h_DoxMinus","H3p3_48h_DoxPlus",
"H3K27ac_UI_DoxMinus","H3K27ac_UI_DoxPlus",
"H3K27ac_0h_DoxMinus","H3K27ac_0h_DoxPlus",
"H3K27ac_24h_DoxMinus","H3K27ac_24h_DoxPlus",
"H3K27ac_48h_DoxMinus","H3K27ac_48h_DoxPlus",
"H3K4me3_UI_DoxMinus","H3K4me3_UI_DoxPlus",
"H3K4me3_0h_DoxMinus","H3K4me3_0h_DoxPlus",
"H3K4me3_24h_DoxMinus","H3K4me3_24h_DoxPlus",
"H3K4me3_48h_DoxMinus","H3K4me3_48h_DoxPlus",
"H3K27me3_UI_DoxMinus","H3K27me3_UI_DoxPlus",
"H3K27me3_0h_DoxMinus","H3K27me3_0h_DoxPlus",
"H3K27me3_24h_DoxMinus","H3K27me3_24h_DoxPlus",
"H3K27me3_48h_DoxMinus","H3K27me3_48h_DoxPlus")
group_H3p3 <- c(
"H3p3_UI_DoxMinus","H3p3_UI_DoxPlus",
"H3p3_0h_DoxMinus","H3p3_0h_DoxPlus",
"H3p3_24h_DoxMinus","H3p3_24h_DoxPlus",
"H3p3_48h_DoxMinus","H3p3_48h_DoxPlus")
group_ATAC <- c(
"ATAC_UI_DoxMinus","ATAC_UI_DoxPlus",
"ATAC_48h_DoxMinus","ATAC_48h_DoxPlus")
samples <- c(
"ATAC_UI_DoxMinus_1","ATAC_UI_DoxMinus_2","ATAC_UI_DoxMinus_3","ATAC_UI_DoxMinus_4",
"ATAC_UI_DoxPlus_1","ATAC_UI_DoxPlus_2","ATAC_UI_DoxPlus_3","ATAC_UI_DoxPlus_4",
"ATAC_48h_DoxMinus_1","ATAC_48h_DoxMinus_2","ATAC_48h_DoxMinus_3","ATAC_48h_DoxMinus_4",
"ATAC_48h_DoxPlus_1","ATAC_48h_DoxPlus_2","ATAC_48h_DoxPlus_3","ATAC_48h_DoxPlus_4",
"H3p3_UI_DoxMinus_1","H3p3_UI_DoxMinus_2","H3p3_UI_DoxPlus_1","H3p3_UI_DoxPlus_2",
"H3p3_0h_DoxMinus_1","H3p3_0h_DoxMinus_2","H3p3_0h_DoxPlus_1","H3p3_0h_DoxPlus_2",
"H3p3_24h_DoxMinus_1","H3p3_24h_DoxMinus_2","H3p3_24h_DoxPlus_1","H3p3_24h_DoxPlus_2",
"H3p3_48h_DoxMinus_1","H3p3_48h_DoxMinus_2","H3p3_48h_DoxPlus_1","H3p3_48h_DoxPlus_2",
"H3K27ac_UI_DoxMinus_1","H3K27ac_UI_DoxMinus_2","H3K27ac_UI_DoxPlus_1","H3K27ac_UI_DoxPlus_2",
"H3K27ac_0h_DoxMinus_1","H3K27ac_0h_DoxMinus_2","H3K27ac_0h_DoxPlus_1","H3K27ac_0h_DoxPlus_2",
"H3K27ac_24h_DoxMinus_1","H3K27ac_24h_DoxMinus_2","H3K27ac_24h_DoxPlus_1","H3K27ac_24h_DoxPlus_2",
"H3K27ac_48h_DoxMinus_1","H3K27ac_48h_DoxMinus_2","H3K27ac_48h_DoxPlus_1","H3K27ac_48h_DoxPlus_2",
"H3K4me3_UI_DoxMinus_1","H3K4me3_UI_DoxMinus_2","H3K4me3_UI_DoxPlus_1","H3K4me3_UI_DoxPlus_2",
"H3K4me3_0h_DoxMinus_1","H3K4me3_0h_DoxMinus_2","H3K4me3_0h_DoxPlus_1","H3K4me3_0h_DoxPlus_2",
"H3K4me3_24h_DoxMinus_1","H3K4me3_24h_DoxMinus_2","H3K4me3_24h_DoxPlus_1","H3K4me3_24h_DoxPlus_2",
"H3K4me3_48h_DoxMinus_1","H3K4me3_48h_DoxMinus_2","H3K4me3_48h_DoxPlus_1","H3K4me3_48h_DoxPlus_2",
"H3K27me3_UI_DoxMinus_1","H3K27me3_UI_DoxMinus_2","H3K27me3_UI_DoxPlus_1","H3K27me3_UI_DoxPlus_2",
"H3K27me3_0h_DoxMinus_1","H3K27me3_0h_DoxMinus_2","H3K27me3_0h_DoxPlus_1","H3K27me3_0h_DoxPlus_2",
"H3K27me3_24h_DoxMinus_1","H3K27me3_24h_DoxMinus_2","H3K27me3_24h_DoxPlus_1","H3K27me3_24h_DoxPlus_2",
"H3K27me3_48h_DoxMinus_1","H3K27me3_48h_DoxMinus_2","H3K27me3_48h_DoxPlus_1","H3K27me3_48h_DoxPlus_2")
samples_ATAC <- c(
"ATAC_UI_DoxMinus_1","ATAC_UI_DoxMinus_2","ATAC_UI_DoxMinus_3","ATAC_UI_DoxMinus_4",
"ATAC_UI_DoxPlus_1","ATAC_UI_DoxPlus_2","ATAC_UI_DoxPlus_3","ATAC_UI_DoxPlus_4",
"ATAC_48h_DoxMinus_1","ATAC_48h_DoxMinus_2","ATAC_48h_DoxMinus_3","ATAC_48h_DoxMinus_4",
"ATAC_48h_DoxPlus_1","ATAC_48h_DoxPlus_2","ATAC_48h_DoxPlus_3","ATAC_48h_DoxPlus_4")
samples_H3p3 <- c(
"H3p3_UI_DoxMinus_1","H3p3_UI_DoxMinus_2","H3p3_UI_DoxPlus_1","H3p3_UI_DoxPlus_2",
"H3p3_0h_DoxMinus_1","H3p3_0h_DoxMinus_2","H3p3_0h_DoxPlus_1","H3p3_0h_DoxPlus_2",
"H3p3_24h_DoxMinus_1","H3p3_24h_DoxMinus_2","H3p3_24h_DoxPlus_1","H3p3_24h_DoxPlus_2",
"H3p3_48h_DoxMinus_1","H3p3_48h_DoxMinus_2","H3p3_48h_DoxPlus_1","H3p3_48h_DoxPlus_2")
samples_H3K27ac <- c(
"H3K27ac_UI_DoxMinus_1","H3K27ac_UI_DoxMinus_2","H3K27ac_UI_DoxPlus_1","H3K27ac_UI_DoxPlus_2",
"H3K27ac_0h_DoxMinus_1","H3K27ac_0h_DoxMinus_2","H3K27ac_0h_DoxPlus_1","H3K27ac_0h_DoxPlus_2",
"H3K27ac_24h_DoxMinus_1","H3K27ac_24h_DoxMinus_2","H3K27ac_24h_DoxPlus_1","H3K27ac_24h_DoxPlus_2",
"H3K27ac_48h_DoxMinus_1","H3K27ac_48h_DoxMinus_2","H3K27ac_48h_DoxPlus_1","H3K27ac_48h_DoxPlus_2")
samples_H3K4me3 <- c(
"H3K4me3_UI_DoxMinus_1","H3K4me3_UI_DoxMinus_2","H3K4me3_UI_DoxPlus_1","H3K4me3_UI_DoxPlus_2",
"H3K4me3_0h_DoxMinus_1","H3K4me3_0h_DoxMinus_2","H3K4me3_0h_DoxPlus_1","H3K4me3_0h_DoxPlus_2",
"H3K4me3_24h_DoxMinus_1","H3K4me3_24h_DoxMinus_2","H3K4me3_24h_DoxPlus_1","H3K4me3_24h_DoxPlus_2",
"H3K4me3_48h_DoxMinus_1","H3K4me3_48h_DoxMinus_2","H3K4me3_48h_DoxPlus_1","H3K4me3_48h_DoxPlus_2")
samples_H3K27me3 <- c(
"H3K27me3_UI_DoxMinus_1","H3K27me3_UI_DoxMinus_2","H3K27me3_UI_DoxPlus_1","H3K27me3_UI_DoxPlus_2",
"H3K27me3_0h_DoxMinus_1","H3K27me3_0h_DoxMinus_2","H3K27me3_0h_DoxPlus_1","H3K27me3_0h_DoxPlus_2",
"H3K27me3_24h_DoxMinus_1","H3K27me3_24h_DoxMinus_2","H3K27me3_24h_DoxPlus_1","H3K27me3_24h_DoxPlus_2",
"H3K27me3_48h_DoxMinus_1","H3K27me3_48h_DoxMinus_2","H3K27me3_48h_DoxPlus_1","H3K27me3_48h_DoxPlus_2")
f_sample <- function(x) x %>% mutate(sample=factor(sample, samples))
f_group <- function(x) x %>% mutate(group=factor(group, groups))
# filter(sample!="Doxminus_D48_ATAC_1") => filter((sample!="Doxminus_D48_ATAC_1")&(rep!="lot4")) (2020 0114修正)
#def_list <- def_list %>% f_sample %>% f_group
####
def_list_select <- def_list
def_list_select_0 <- def_list %>% dplyr::select(-"file",-"multicov_No") %>% filter(seq=="ATAC")
def_list_select_1 <- def_list %>% dplyr::select(-"file",-"multicov_No") %>% filter(seq=="H3p3")
def_list_select_2 <- def_list %>% dplyr::select(-"file",-"multicov_No") %>% filter(seq=="H3K27ac")
def_list_select_3 <- def_list %>% dplyr::select(-"file",-"multicov_No") %>% filter(seq=="H3K4me3")
def_list_select_4 <- def_list %>% dplyr::select(-"file",-"multicov_No") %>% filter(seq=="H3K27me3")
#def_list_select <- def_list %>% filter(!!use) #使わないサンプルを削除(20190917)
#%>%
#mutate(time1 = time, rep1=rep) %>% unite(time1,rep1,col="time_replicate") %>% #mutate(time_replicate=factor(time_replicate,c("UI_1","UI_2","D1","D48_lot2")))
# narrowpeak の場合 (型を指定) 200612modif 200615modif
#narrow_colnames <- c("chr","start","end","name","score","strand","singnalValue","pValue","qValue","peak","chr0","start0","end0")
merge_colnames <- c("chr","TSSstart","TSSend","ens_gene","score","strand","TSS","Start","End")
matome0 <- readr::read_tsv(filepath_summary, col_names = c(merge_colnames, def_list$sample)) %>% mutate_if(names(.) %in% c("start","end","TSS","Start","End",def_list$sample), as.integer)
Parsed with column specification:
cols(
.default = col_double(),
chr = col_character(),
ens_gene = col_character(),
score = col_character(),
strand = col_character()
)
See spec(...) for full column specifications.
matome0 %>% dplyr::select("chr","TSSstart","TSSend","ens_gene") %>% unique() # # A tibble: 105,166 x 4
#matome0 %>% group_by(chr,ens_gene,strand,Start,End) %>% unique() # # A tibble: 105,166 x 4
#matome <- matome0 %>% group_by(chr,start,end,Name) %>% dplyr::top_n(1, singnalValue) %>% dplyr::top_n(1, peak) %>% dplyr::ungroup()
#matome <- matome0 %>% filter_at(def_list$sample,any_vars(. > 0))
matome1 <- matome0 %>% gather("sample", "count", -("chr":"End")) %>% filter(sample %in% def_list_select$sample) %>% left_join(def_list_select, .,by = "sample")
# %>% dplyr::select(-"score",-"strand",-"singnalValue",-"pValue",-"qValue",-"peak")
#---- 確認 ----#
matome0 %>% nrow()
[1] 133122
filename <- gsub(".txt","__count.csv",basename(filepath_summary)) #geneにつき複数領域
matome0 %>% readr::write_csv(filename)
print(filename)
[1] "ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__TSS_pm5kb_20200624.count__count.csv"
#filename <- gsub(".txt","__select_nameonly.csv",basename(filepath_def))
#matome %>% dplyr::select("chr","start","end","name") %>% readr::write_csv(filename)
#print(filename)
#matome %>% dplyr::select("chr","start","end","name","score","Name","strand","singnalValue","pValue","qValue","peak","chr0","start0","end0") %>% readr::write_csv(filename)
#print(filename)
#matome00 <- matome0 %>% group_by(chr,start,end,Name) %>% summarise(max_score=max(score), paste(Name, collapse="")) #dplyr::top_n(score,1)
#matome <- matome0 %>% filter(name %in% matome00$name)
#matome %>% mutate(name1=name) %>% filter(grepl("[a-u]$",name)) %>% mutate(Name=gsub("[a-u]$","",name))
#mat_select <- matome %>% dplyr::select(chr,start,end,name,Name)
#annotate_bed <- partial(right_join,mat_select,by="ens_gene") #2gunで使う
#matome1_number <- matome1 %>% mutate(position = row_number())
#matome1_plus <- matome1_number %>% filter(strand=="+")
#matome1_minus <- matome1_number %>% filter(strand=="-")
matome0_number <- matome0 %>% mutate(position = row_number())
nrow(matome0_number)
[1] 133122
matome0_plus <- matome0_number %>% filter(strand=="+")
matome0_minus <- matome0_number %>% filter(strand=="-")
matome0_plus_o <- matome0_plus %>% group_by(chr,ens_gene) %>% dplyr::top_n(-1,TSS) #低
matome0_minus_o <- matome0_minus %>% group_by(chr,ens_gene) %>% dplyr::top_n(1,TSS) #高
matome0_o <- dplyr::bind_rows(matome0_plus_o, matome0_minus_o) %>% arrange(position)
nrow(matome0_o)
[1] 55487
##----- 確認 ---------##
colnames(matome0_o)
[1] "chr" "TSSstart" "TSSend" "ens_gene" "score"
[6] "strand" "TSS" "Start" "End" "H3p3_UI_DoxMinus_1"
[11] "H3p3_UI_DoxMinus_2" "H3p3_UI_DoxPlus_1" "H3p3_UI_DoxPlus_2" "H3p3_0h_DoxMinus_1" "H3p3_0h_DoxMinus_2"
[16] "H3p3_0h_DoxPlus_1" "H3p3_0h_DoxPlus_2" "H3p3_24h_DoxMinus_1" "H3p3_24h_DoxMinus_2" "H3p3_24h_DoxPlus_1"
[21] "H3p3_24h_DoxPlus_2" "H3p3_48h_DoxMinus_1" "H3p3_48h_DoxMinus_2" "H3p3_48h_DoxPlus_1" "H3p3_48h_DoxPlus_2"
[26] "H3K27ac_UI_DoxMinus_1" "H3K27ac_UI_DoxMinus_2" "H3K27ac_UI_DoxPlus_1" "H3K27ac_UI_DoxPlus_2" "H3K27ac_0h_DoxMinus_1"
[31] "H3K27ac_0h_DoxMinus_2" "H3K27ac_0h_DoxPlus_1" "H3K27ac_0h_DoxPlus_2" "H3K27ac_24h_DoxMinus_1" "H3K27ac_24h_DoxMinus_2"
[36] "H3K27ac_24h_DoxPlus_1" "H3K27ac_24h_DoxPlus_2" "H3K27ac_48h_DoxMinus_1" "H3K27ac_48h_DoxMinus_2" "H3K27ac_48h_DoxPlus_1"
[41] "H3K27ac_48h_DoxPlus_2" "H3K4me3_UI_DoxMinus_1" "H3K4me3_UI_DoxMinus_2" "H3K4me3_UI_DoxPlus_1" "H3K4me3_UI_DoxPlus_2"
[46] "H3K4me3_0h_DoxMinus_1" "H3K4me3_0h_DoxMinus_2" "H3K4me3_0h_DoxPlus_1" "H3K4me3_0h_DoxPlus_2" "H3K4me3_24h_DoxMinus_1"
[51] "H3K4me3_24h_DoxMinus_2" "H3K4me3_24h_DoxPlus_1" "H3K4me3_24h_DoxPlus_2" "H3K4me3_48h_DoxMinus_1" "H3K4me3_48h_DoxMinus_2"
[56] "H3K4me3_48h_DoxPlus_1" "H3K4me3_48h_DoxPlus_2" "H3K27me3_UI_DoxMinus_1" "H3K27me3_UI_DoxMinus_2" "H3K27me3_UI_DoxPlus_1"
[61] "H3K27me3_UI_DoxPlus_2" "H3K27me3_0h_DoxMinus_1" "H3K27me3_0h_DoxMinus_2" "H3K27me3_0h_DoxPlus_1" "H3K27me3_0h_DoxPlus_2"
[66] "H3K27me3_24h_DoxMinus_1" "H3K27me3_24h_DoxMinus_2" "H3K27me3_24h_DoxPlus_1" "H3K27me3_24h_DoxPlus_2" "H3K27me3_48h_DoxMinus_1"
[71] "H3K27me3_48h_DoxMinus_2" "H3K27me3_48h_DoxPlus_1" "H3K27me3_48h_DoxPlus_2" "ATAC_UI_DoxMinus_1" "ATAC_UI_DoxMinus_2"
[76] "ATAC_UI_DoxMinus_3" "ATAC_UI_DoxMinus_4" "ATAC_UI_DoxPlus_1" "ATAC_UI_DoxPlus_2" "ATAC_UI_DoxPlus_3"
[81] "ATAC_UI_DoxPlus_4" "ATAC_48h_DoxMinus_1" "ATAC_48h_DoxMinus_2" "ATAC_48h_DoxMinus_3" "ATAC_48h_DoxMinus_4"
[86] "ATAC_48h_DoxPlus_1" "ATAC_48h_DoxPlus_2" "ATAC_48h_DoxPlus_3" "ATAC_48h_DoxPlus_4" "position"
nrow(matome0_o)
[1] 55487
filename <- gsub(".txt","__count_firstTSS.csv",basename(filepath_summary)) #geneにつき1領域
matome0_o %>% readr::write_csv(filename)
print(filename)
[1] "ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__TSS_pm5kb_20200624.count__count_firstTSS.csv"
##--------------------##
#matome0_plus %>% filter(ens_gene =="ENSMUSG00000000037") %>% dplyr::select(TSSstart,position)
#matome0_plus_o %>% filter(ens_gene =="ENSMUSG00000000037") %>% dplyr::select(TSSstart,position)
#%>% group_by(chr,TSSstart,TSSend,ens_gene,score,strand,TSS,Start,End,seq)
#matome1_number %>% filter(!(strand=="+"|strand=="-"))
matome0_s <- matome0_o %>% dplyr::select(chr,ens_gene,TSSstart,TSSend,score,strand,TSS,Start,End,position,all_of(samples)) %>% filter(chr!="chrM")
nrow(matome0_s)
[1] 55450
#matome0_s <- matome0_o %>% filter(ens_gene %in% FC_rank_all_BRBlist$ens_gene) %>% dplyr::select(chr,ens_gene,TSSstart,TSSend,score,strand,TSS,Start,End,position,all_of(samples))
nrow(matome0_s)
[1] 55450
#matome0_s1 <- matome0_s %>% left_join(FC_rank_all_BRBlist %>% dplyr::select(ens_gene,log2FoldChange,Rank))
matome5 <- matome0_s %>% dplyr::select(chr,ens_gene,position,all_of(samples)) %>% ungroup()
##----- 確認 ---------##
colnames(matome0_s)
[1] "chr" "ens_gene" "TSSstart" "TSSend" "score"
[6] "strand" "TSS" "Start" "End" "position"
[11] "ATAC_UI_DoxMinus_1" "ATAC_UI_DoxMinus_2" "ATAC_UI_DoxMinus_3" "ATAC_UI_DoxMinus_4" "ATAC_UI_DoxPlus_1"
[16] "ATAC_UI_DoxPlus_2" "ATAC_UI_DoxPlus_3" "ATAC_UI_DoxPlus_4" "ATAC_48h_DoxMinus_1" "ATAC_48h_DoxMinus_2"
[21] "ATAC_48h_DoxMinus_3" "ATAC_48h_DoxMinus_4" "ATAC_48h_DoxPlus_1" "ATAC_48h_DoxPlus_2" "ATAC_48h_DoxPlus_3"
[26] "ATAC_48h_DoxPlus_4" "H3p3_UI_DoxMinus_1" "H3p3_UI_DoxMinus_2" "H3p3_UI_DoxPlus_1" "H3p3_UI_DoxPlus_2"
[31] "H3p3_0h_DoxMinus_1" "H3p3_0h_DoxMinus_2" "H3p3_0h_DoxPlus_1" "H3p3_0h_DoxPlus_2" "H3p3_24h_DoxMinus_1"
[36] "H3p3_24h_DoxMinus_2" "H3p3_24h_DoxPlus_1" "H3p3_24h_DoxPlus_2" "H3p3_48h_DoxMinus_1" "H3p3_48h_DoxMinus_2"
[41] "H3p3_48h_DoxPlus_1" "H3p3_48h_DoxPlus_2" "H3K27ac_UI_DoxMinus_1" "H3K27ac_UI_DoxMinus_2" "H3K27ac_UI_DoxPlus_1"
[46] "H3K27ac_UI_DoxPlus_2" "H3K27ac_0h_DoxMinus_1" "H3K27ac_0h_DoxMinus_2" "H3K27ac_0h_DoxPlus_1" "H3K27ac_0h_DoxPlus_2"
[51] "H3K27ac_24h_DoxMinus_1" "H3K27ac_24h_DoxMinus_2" "H3K27ac_24h_DoxPlus_1" "H3K27ac_24h_DoxPlus_2" "H3K27ac_48h_DoxMinus_1"
[56] "H3K27ac_48h_DoxMinus_2" "H3K27ac_48h_DoxPlus_1" "H3K27ac_48h_DoxPlus_2" "H3K4me3_UI_DoxMinus_1" "H3K4me3_UI_DoxMinus_2"
[61] "H3K4me3_UI_DoxPlus_1" "H3K4me3_UI_DoxPlus_2" "H3K4me3_0h_DoxMinus_1" "H3K4me3_0h_DoxMinus_2" "H3K4me3_0h_DoxPlus_1"
[66] "H3K4me3_0h_DoxPlus_2" "H3K4me3_24h_DoxMinus_1" "H3K4me3_24h_DoxMinus_2" "H3K4me3_24h_DoxPlus_1" "H3K4me3_24h_DoxPlus_2"
[71] "H3K4me3_48h_DoxMinus_1" "H3K4me3_48h_DoxMinus_2" "H3K4me3_48h_DoxPlus_1" "H3K4me3_48h_DoxPlus_2" "H3K27me3_UI_DoxMinus_1"
[76] "H3K27me3_UI_DoxMinus_2" "H3K27me3_UI_DoxPlus_1" "H3K27me3_UI_DoxPlus_2" "H3K27me3_0h_DoxMinus_1" "H3K27me3_0h_DoxMinus_2"
[81] "H3K27me3_0h_DoxPlus_1" "H3K27me3_0h_DoxPlus_2" "H3K27me3_24h_DoxMinus_1" "H3K27me3_24h_DoxMinus_2" "H3K27me3_24h_DoxPlus_1"
[86] "H3K27me3_24h_DoxPlus_2" "H3K27me3_48h_DoxMinus_1" "H3K27me3_48h_DoxMinus_2" "H3K27me3_48h_DoxPlus_1" "H3K27me3_48h_DoxPlus_2"
nrow(matome0_s)
[1] 55450
filename <- gsub(".txt","__count_firstTSS_select.csv",basename(filepath_summary)) #geneにつき1領域、かつBRBでnormalized countがあるもの
matome0_s %>% readr::write_csv(filename)
print(filename)
[1] "ChIL01100111_ATAC0049L1__3T3_EGFP18_Dox__TSS_pm5kb_20200624.count__count_firstTSS_select.csv"
##--------------------##
annotate_TSS <- partial(right_join,dplyr::select(matome0_s,chr,ens_gene,TSSstart,TSSend,score,strand,TSS,Start,End,position),by="ens_gene") #2gunで使う
separate matrix
nrow(matome5)
[1] 55450
colnames(matome5)
[1] "chr" "ens_gene" "position" "ATAC_UI_DoxMinus_1" "ATAC_UI_DoxMinus_2"
[6] "ATAC_UI_DoxMinus_3" "ATAC_UI_DoxMinus_4" "ATAC_UI_DoxPlus_1" "ATAC_UI_DoxPlus_2" "ATAC_UI_DoxPlus_3"
[11] "ATAC_UI_DoxPlus_4" "ATAC_48h_DoxMinus_1" "ATAC_48h_DoxMinus_2" "ATAC_48h_DoxMinus_3" "ATAC_48h_DoxMinus_4"
[16] "ATAC_48h_DoxPlus_1" "ATAC_48h_DoxPlus_2" "ATAC_48h_DoxPlus_3" "ATAC_48h_DoxPlus_4" "H3p3_UI_DoxMinus_1"
[21] "H3p3_UI_DoxMinus_2" "H3p3_UI_DoxPlus_1" "H3p3_UI_DoxPlus_2" "H3p3_0h_DoxMinus_1" "H3p3_0h_DoxMinus_2"
[26] "H3p3_0h_DoxPlus_1" "H3p3_0h_DoxPlus_2" "H3p3_24h_DoxMinus_1" "H3p3_24h_DoxMinus_2" "H3p3_24h_DoxPlus_1"
[31] "H3p3_24h_DoxPlus_2" "H3p3_48h_DoxMinus_1" "H3p3_48h_DoxMinus_2" "H3p3_48h_DoxPlus_1" "H3p3_48h_DoxPlus_2"
[36] "H3K27ac_UI_DoxMinus_1" "H3K27ac_UI_DoxMinus_2" "H3K27ac_UI_DoxPlus_1" "H3K27ac_UI_DoxPlus_2" "H3K27ac_0h_DoxMinus_1"
[41] "H3K27ac_0h_DoxMinus_2" "H3K27ac_0h_DoxPlus_1" "H3K27ac_0h_DoxPlus_2" "H3K27ac_24h_DoxMinus_1" "H3K27ac_24h_DoxMinus_2"
[46] "H3K27ac_24h_DoxPlus_1" "H3K27ac_24h_DoxPlus_2" "H3K27ac_48h_DoxMinus_1" "H3K27ac_48h_DoxMinus_2" "H3K27ac_48h_DoxPlus_1"
[51] "H3K27ac_48h_DoxPlus_2" "H3K4me3_UI_DoxMinus_1" "H3K4me3_UI_DoxMinus_2" "H3K4me3_UI_DoxPlus_1" "H3K4me3_UI_DoxPlus_2"
[56] "H3K4me3_0h_DoxMinus_1" "H3K4me3_0h_DoxMinus_2" "H3K4me3_0h_DoxPlus_1" "H3K4me3_0h_DoxPlus_2" "H3K4me3_24h_DoxMinus_1"
[61] "H3K4me3_24h_DoxMinus_2" "H3K4me3_24h_DoxPlus_1" "H3K4me3_24h_DoxPlus_2" "H3K4me3_48h_DoxMinus_1" "H3K4me3_48h_DoxMinus_2"
[66] "H3K4me3_48h_DoxPlus_1" "H3K4me3_48h_DoxPlus_2" "H3K27me3_UI_DoxMinus_1" "H3K27me3_UI_DoxMinus_2" "H3K27me3_UI_DoxPlus_1"
[71] "H3K27me3_UI_DoxPlus_2" "H3K27me3_0h_DoxMinus_1" "H3K27me3_0h_DoxMinus_2" "H3K27me3_0h_DoxPlus_1" "H3K27me3_0h_DoxPlus_2"
[76] "H3K27me3_24h_DoxMinus_1" "H3K27me3_24h_DoxMinus_2" "H3K27me3_24h_DoxPlus_1" "H3K27me3_24h_DoxPlus_2" "H3K27me3_48h_DoxMinus_1"
[81] "H3K27me3_48h_DoxMinus_2" "H3K27me3_48h_DoxPlus_1" "H3K27me3_48h_DoxPlus_2"
X <- matome5 %>% dplyr::select(all_of(samples)) %>% as.matrix
rownames(X) <- matome5$ens_gene
###--- DESeq2によりnormalized count (必要なサンプルのみ) ---###
model<- ~group
dds_ATAC <- DESeq2::DESeqDataSetFromMatrix(X[,def_list_select_0$sample],def_list_select_0,model) #ATAC
some variables in design formula are characters, converting to factors
dds_H3p3 <- DESeq2::DESeqDataSetFromMatrix(X[,def_list_select_1$sample],def_list_select_1,model) #H3p3
some variables in design formula are characters, converting to factors
dds_H3K27ac <- DESeq2::DESeqDataSetFromMatrix(X[,def_list_select_2$sample],def_list_select_2,model) #H3K27ac
some variables in design formula are characters, converting to factors
dds_H3K4me3 <- DESeq2::DESeqDataSetFromMatrix(X[,def_list_select_3$sample],def_list_select_3,model) #H3K4me3
some variables in design formula are characters, converting to factors
dds_H3K27me3 <- DESeq2::DESeqDataSetFromMatrix(X[,def_list_select_4$sample],def_list_select_4,model) #H3K27me3
some variables in design formula are characters, converting to factors
#model_2gun <- ~group
#dds_2gun <- DESeq2::DESeqDataSetFromMatrix(X[,def_list_select$sample],def_list_select,model_2gun)
dds_ATAC <- DESeq2::DESeq(dds_ATAC)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
dds_H3p3 <- DESeq2::DESeq(dds_H3p3)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
dds_H3K27ac <- DESeq2::DESeq(dds_H3K27ac)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing
dds_H3K4me3 <- DESeq2::DESeq(dds_H3K4me3)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing
dds_H3K27me3 <- DESeq2::DESeq(dds_H3K27me3)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing
#keep <- rowSums(counts(dds)) >= 10 #low countは削る方法
#dds <- dds[keep,] #low countは削る方法
DESeq2::sizeFactors(dds_ATAC) %>%
{tibble(sample=names(.),sizeFactor=.)} %>%
ggplot(aes(sample,sizeFactor)) + theme_minimal() +
geom_bar(stat="identity") + coord_flip()
DESeq2::plotDispEsts(dds_ATAC)
DESeq2::sizeFactors(dds_H3p3) %>%
{tibble(sample=names(.),sizeFactor=.)} %>%
ggplot(aes(sample,sizeFactor)) + theme_minimal() +
geom_bar(stat="identity") + coord_flip()
DESeq2::plotDispEsts(dds_H3p3)
DESeq2::sizeFactors(dds_H3K27ac) %>%
{tibble(sample=names(.),sizeFactor=.)} %>%
ggplot(aes(sample,sizeFactor)) + theme_minimal() +
geom_bar(stat="identity") + coord_flip()
DESeq2::plotDispEsts(dds_H3K27ac)
DESeq2::sizeFactors(dds_H3K4me3) %>%
{tibble(sample=names(.),sizeFactor=.)} %>%
ggplot(aes(sample,sizeFactor)) + theme_minimal() +
geom_bar(stat="identity") + coord_flip()
DESeq2::plotDispEsts(dds_H3K4me3)
DESeq2::sizeFactors(dds_H3K27me3) %>%
{tibble(sample=names(.),sizeFactor=.)} %>%
ggplot(aes(sample,sizeFactor)) + theme_minimal() +
geom_bar(stat="identity") + coord_flip()
DESeq2::plotDispEsts(dds_H3K27me3)
normalized count listを書き出し
# ATAC
dds_ATAC <- DESeq2::estimateSizeFactors(dds_ATAC)
norm_ATAC <- DESeq2::counts(dds_ATAC,normalized=TRUE) #DEGを取った後のクラスタリングに使う。
normalizedcount_ATAC <- as.data.frame(norm_ATAC) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_0$sample))
filename <- paste(csvfilepath,"_normcount_ATAC.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_ATAC.csv"
readr::write_csv(normalizedcount_ATAC, filename)
nrow(normalizedcount_ATAC)
[1] 55450
ncol(normalizedcount_ATAC)
[1] 17
norm_gene_ATAC <- normalizedcount_ATAC %>% inner_join(ensemble) %>% dplyr::select("ens_gene","ext_gene", "biotype","chr", all_of(def_list_select_0$sample))
Joining, by = "ens_gene"
filename <- paste(csvfilepath,"_normcount_ATAC_genename.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_ATAC_genename.csv"
readr::write_csv(norm_gene_ATAC, filename)
nrow(norm_gene_ATAC)
[1] 55197
ncol(norm_gene_ATAC)
[1] 20
# H3p3
dds_H3p3 <- DESeq2::estimateSizeFactors(dds_H3p3)
norm_H3p3 <- DESeq2::counts(dds_H3p3,normalized=TRUE) #DEGを取った後のクラスタリングに使う。
normalizedcount_H3p3 <- as.data.frame(norm_H3p3) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_1$sample))
filename <- paste(csvfilepath,"_normcount_H3p3.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3p3.csv"
readr::write_csv(normalizedcount_H3p3, filename)
nrow(normalizedcount_H3p3)
[1] 55450
ncol(normalizedcount_H3p3)
[1] 17
norm_gene_H3p3 <- normalizedcount_H3p3 %>% inner_join(ensemble) %>% dplyr::select("ens_gene","ext_gene", "biotype","chr", all_of(def_list_select_1$sample))
Joining, by = "ens_gene"
filename <- paste(csvfilepath,"_normcount_H3p3_genename.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3p3_genename.csv"
readr::write_csv(norm_gene_H3p3, filename)
nrow(norm_gene_H3p3)
[1] 55197
ncol(norm_gene_H3p3)
[1] 20
# H3K27ac
dds_H3K27ac <- DESeq2::estimateSizeFactors(dds_H3K27ac)
norm_H3K27ac <- DESeq2::counts(dds_H3K27ac,normalized=TRUE) #DEGを取った後のクラスタリングに使う。
normalizedcount_H3K27ac <- as.data.frame(norm_H3K27ac) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_2$sample))
filename <- paste(csvfilepath,"_normcount_H3K27ac.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3K27ac.csv"
readr::write_csv(normalizedcount_H3K27ac, filename)
nrow(normalizedcount_H3K27ac)
[1] 55450
ncol(normalizedcount_H3K27ac)
[1] 17
norm_gene_H3K27ac <- normalizedcount_H3K27ac %>% inner_join(ensemble) %>% dplyr::select("ens_gene","ext_gene", "biotype","chr", all_of(def_list_select_2$sample))
Joining, by = "ens_gene"
filename <- paste(csvfilepath,"_normcount_H3K27ac_genename.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3K27ac_genename.csv"
readr::write_csv(norm_gene_H3K27ac, filename)
nrow(norm_gene_H3K27ac)
[1] 55197
ncol(norm_gene_H3K27ac)
[1] 20
# H3K4me3
dds_H3K4me3 <- DESeq2::estimateSizeFactors(dds_H3K4me3)
norm_H3K4me3 <- DESeq2::counts(dds_H3K4me3,normalized=TRUE) #DEGを取った後のクラスタリングに使う。
normalizedcount_H3K4me3 <- as.data.frame(norm_H3K4me3) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_3$sample))
filename <- paste(csvfilepath,"_normcount_H3K4me3.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3K4me3.csv"
readr::write_csv(normalizedcount_H3K4me3, filename)
nrow(normalizedcount_H3K4me3)
[1] 55450
ncol(normalizedcount_H3K4me3)
[1] 17
norm_gene_H3K4me3 <- normalizedcount_H3K4me3 %>% inner_join(ensemble) %>% dplyr::select("ens_gene","ext_gene", "biotype","chr", all_of(def_list_select_3$sample))
Joining, by = "ens_gene"
filename <- paste(csvfilepath,"_normcount_H3K4me3_genename.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3K4me3_genename.csv"
readr::write_csv(norm_gene_H3K4me3, filename)
nrow(norm_gene_H3K4me3)
[1] 55197
ncol(norm_gene_H3K4me3)
[1] 20
# H3K27me3
dds_H3K27me3 <- DESeq2::estimateSizeFactors(dds_H3K27me3)
norm_H3K27me3 <- DESeq2::counts(dds_H3K27me3,normalized=TRUE) #DEGを取った後のクラスタリングに使う。
normalizedcount_H3K27me3 <- as.data.frame(norm_H3K27me3) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_4$sample))
filename <- paste(csvfilepath,"_normcount_H3K27me3.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3K27me3.csv"
readr::write_csv(normalizedcount_H3K27me3, filename)
nrow(normalizedcount_H3K27me3)
[1] 55450
ncol(normalizedcount_H3K27me3)
[1] 17
norm_gene_H3K27me3 <- normalizedcount_H3K27me3 %>% inner_join(ensemble) %>% dplyr::select("ens_gene","ext_gene", "biotype","chr", all_of(def_list_select_4$sample))
Joining, by = "ens_gene"
filename <- paste(csvfilepath,"_normcount_H3K27me3_genename.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_H3K27me3_genename.csv"
readr::write_csv(norm_gene_H3K27me3, filename)
nrow(norm_gene_H3K27me3)
[1] 55197
ncol(norm_gene_H3K27me3)
[1] 20
# bind norm count
normalizedcount <- normalizedcount_ATAC %>% inner_join(normalizedcount_H3p3) %>% inner_join(normalizedcount_H3K27ac) %>% inner_join(normalizedcount_H3K4me3) %>% inner_join(normalizedcount_H3K27me3)
Joining, by = "ens_gene"
Joining, by = "ens_gene"
Joining, by = "ens_gene"
Joining, by = "ens_gene"
norm_gene <- norm_gene_ATAC %>% inner_join(norm_gene_H3p3) %>% inner_join(norm_gene_H3K27ac) %>% inner_join(norm_gene_H3K4me3) %>% inner_join(norm_gene_H3K27me3)
Joining, by = c("ens_gene", "ext_gene", "biotype", "chr")
Joining, by = c("ens_gene", "ext_gene", "biotype", "chr")
Joining, by = c("ens_gene", "ext_gene", "biotype", "chr")
Joining, by = c("ens_gene", "ext_gene", "biotype", "chr")
print(norm_gene)
colnames(norm_gene)
[1] "ens_gene" "ext_gene" "biotype" "chr" "ATAC_UI_DoxMinus_1"
[6] "ATAC_UI_DoxMinus_2" "ATAC_UI_DoxMinus_3" "ATAC_UI_DoxMinus_4" "ATAC_UI_DoxPlus_1" "ATAC_UI_DoxPlus_2"
[11] "ATAC_UI_DoxPlus_3" "ATAC_UI_DoxPlus_4" "ATAC_48h_DoxMinus_1" "ATAC_48h_DoxMinus_2" "ATAC_48h_DoxMinus_3"
[16] "ATAC_48h_DoxMinus_4" "ATAC_48h_DoxPlus_1" "ATAC_48h_DoxPlus_2" "ATAC_48h_DoxPlus_3" "ATAC_48h_DoxPlus_4"
[21] "H3p3_UI_DoxMinus_1" "H3p3_UI_DoxMinus_2" "H3p3_UI_DoxPlus_1" "H3p3_UI_DoxPlus_2" "H3p3_0h_DoxMinus_1"
[26] "H3p3_0h_DoxMinus_2" "H3p3_0h_DoxPlus_1" "H3p3_0h_DoxPlus_2" "H3p3_24h_DoxMinus_1" "H3p3_24h_DoxMinus_2"
[31] "H3p3_24h_DoxPlus_1" "H3p3_24h_DoxPlus_2" "H3p3_48h_DoxMinus_1" "H3p3_48h_DoxMinus_2" "H3p3_48h_DoxPlus_1"
[36] "H3p3_48h_DoxPlus_2" "H3K27ac_UI_DoxMinus_1" "H3K27ac_UI_DoxMinus_2" "H3K27ac_UI_DoxPlus_1" "H3K27ac_UI_DoxPlus_2"
[41] "H3K27ac_0h_DoxMinus_1" "H3K27ac_0h_DoxMinus_2" "H3K27ac_0h_DoxPlus_1" "H3K27ac_0h_DoxPlus_2" "H3K27ac_24h_DoxMinus_1"
[46] "H3K27ac_24h_DoxMinus_2" "H3K27ac_24h_DoxPlus_1" "H3K27ac_24h_DoxPlus_2" "H3K27ac_48h_DoxMinus_1" "H3K27ac_48h_DoxMinus_2"
[51] "H3K27ac_48h_DoxPlus_1" "H3K27ac_48h_DoxPlus_2" "H3K4me3_UI_DoxMinus_1" "H3K4me3_UI_DoxMinus_2" "H3K4me3_UI_DoxPlus_1"
[56] "H3K4me3_UI_DoxPlus_2" "H3K4me3_0h_DoxMinus_1" "H3K4me3_0h_DoxMinus_2" "H3K4me3_0h_DoxPlus_1" "H3K4me3_0h_DoxPlus_2"
[61] "H3K4me3_24h_DoxMinus_1" "H3K4me3_24h_DoxMinus_2" "H3K4me3_24h_DoxPlus_1" "H3K4me3_24h_DoxPlus_2" "H3K4me3_48h_DoxMinus_1"
[66] "H3K4me3_48h_DoxMinus_2" "H3K4me3_48h_DoxPlus_1" "H3K4me3_48h_DoxPlus_2" "H3K27me3_UI_DoxMinus_1" "H3K27me3_UI_DoxMinus_2"
[71] "H3K27me3_UI_DoxPlus_1" "H3K27me3_UI_DoxPlus_2" "H3K27me3_0h_DoxMinus_1" "H3K27me3_0h_DoxMinus_2" "H3K27me3_0h_DoxPlus_1"
[76] "H3K27me3_0h_DoxPlus_2" "H3K27me3_24h_DoxMinus_1" "H3K27me3_24h_DoxMinus_2" "H3K27me3_24h_DoxPlus_1" "H3K27me3_24h_DoxPlus_2"
[81] "H3K27me3_48h_DoxMinus_1" "H3K27me3_48h_DoxMinus_2" "H3K27me3_48h_DoxPlus_1" "H3K27me3_48h_DoxPlus_2"
filename <- paste(csvfilepath,"_normcount.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount.csv"
readr::write_csv(normalizedcount, filename)
filename <- paste(csvfilepath,"_normcount_genename.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__normcount_genename.csv"
readr::write_csv(norm_gene, filename)
size factors を書き出し
filename <- paste(csvfilepath,"_sizefactors_ATAC.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__sizefactors_ATAC.csv"
as.data.frame(DESeq2::sizeFactors(dds_ATAC)) %>% tibble::rownames_to_column("sample") %>% readr::write_csv(filename)
filename <- paste(csvfilepath,"_sizefactors_H3p3.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__sizefactors_H3p3.csv"
as.data.frame(DESeq2::sizeFactors(dds_H3p3)) %>% tibble::rownames_to_column("sample") %>% readr::write_csv(filename)
filename <- paste(csvfilepath,"_sizefactors_H3K27ac.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__sizefactors_H3K27ac.csv"
as.data.frame(DESeq2::sizeFactors(dds_H3K27ac)) %>% tibble::rownames_to_column("sample") %>% readr::write_csv(filename)
filename <- paste(csvfilepath,"_sizefactors_H3K4me3.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__sizefactors_H3K4me3.csv"
as.data.frame(DESeq2::sizeFactors(dds_H3K4me3)) %>% tibble::rownames_to_column("sample") %>% readr::write_csv(filename)
filename <- paste(csvfilepath,"_sizefactors_H3K27me3.csv",sep="_")
print(filename)
[1] "TSS_pm5kb_20200624__sizefactors_H3K27me3.csv"
as.data.frame(DESeq2::sizeFactors(dds_H3K27me3)) %>% tibble::rownames_to_column("sample") %>% readr::write_csv(filename)
vst => z score
vsd_ATAC <- DESeq2::vst(dds_ATAC) #normalized countが入っている。(vstかrlog)
Xd_ATAC <- SummarizedExperiment::assay(vsd_ATAC) # 全て選択(200326) 20190920を元に (191024)
Xs_ATAC <- Xd_ATAC %>% t %>% scale %>% t
vsd_H3p3 <- DESeq2::vst(dds_H3p3) #normalized countが入っている。(vstかrlog)
Xd_H3p3 <- SummarizedExperiment::assay(vsd_H3p3) # 全て選択(200326) 20190920を元に (191024)
Xs_H3p3 <- Xd_H3p3 %>% t %>% scale %>% t
vsd_H3K27ac <- DESeq2::vst(dds_H3K27ac) #normalized countが入っている。(vstかrlog)
Xd_H3K27ac <- SummarizedExperiment::assay(vsd_H3K27ac) # 全て選択(200326) 20190920を元に (191024)
Xs_H3K27ac <- Xd_H3K27ac %>% t %>% scale %>% t
vsd_H3K4me3 <- DESeq2::vst(dds_H3K4me3) #normalized countが入っている。(vstかrlog)
Xd_H3K4me3 <- SummarizedExperiment::assay(vsd_H3K4me3) # 全て選択(200326) 20190920を元に (191024)
Xs_H3K4me3 <- Xd_H3K4me3 %>% t %>% scale %>% t
vsd_H3K27me3 <- DESeq2::vst(dds_H3K27me3) #normalized countが入っている。(vstかrlog)
Xd_H3K27me3 <- SummarizedExperiment::assay(vsd_H3K27me3) # 全て選択(200326) 20190920を元に (191024)
Xs_H3K27me3 <- Xd_H3K27me3 %>% t %>% scale %>% t
vsdtrans_ATAC <- as.data.frame(Xd_ATAC) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_0$sample))
filename <- paste(csvfilepath,"_vstrans_ATAC.csv",sep="_")
readr::write_csv(vsdtrans_ATAC, filename)
nrow(vsdtrans_ATAC)
[1] 55450
ncol(vsdtrans_ATAC)
[1] 17
vsdtrans_H3p3 <- as.data.frame(Xd_H3p3) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_1$sample))
filename <- paste(csvfilepath,"_vstrans_H3p3.csv",sep="_")
readr::write_csv(vsdtrans_H3p3, filename)
nrow(vsdtrans_H3p3)
[1] 55450
ncol(vsdtrans_H3p3)
[1] 17
vsdtrans_H3K27ac <- as.data.frame(Xd_H3K27ac) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_2$sample))
filename <- paste(csvfilepath,"_vstrans_H3K27ac.csv",sep="_")
readr::write_csv(vsdtrans_H3K27ac, filename)
nrow(vsdtrans_H3K27ac)
[1] 55450
ncol(vsdtrans_H3K27ac)
[1] 17
vsdtrans_H3K4me3 <- as.data.frame(Xd_H3K4me3) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_3$sample))
filename <- paste(csvfilepath,"_vstrans_H3K4me3.csv",sep="_")
readr::write_csv(vsdtrans_H3K4me3, filename)
nrow(vsdtrans_H3K4me3)
[1] 55450
ncol(vsdtrans_H3K4me3)
[1] 17
vsdtrans_H3K27me3 <- as.data.frame(Xd_H3K27me3) %>% tibble::rownames_to_column("ens_gene") %>% as_tibble %>% dplyr::select("ens_gene", all_of(def_list_select_4$sample))
filename <- paste(csvfilepath,"_vstrans_H3K27me3.csv",sep="_")
readr::write_csv(vsdtrans_H3K27me3, filename)
nrow(vsdtrans_H3K27me3)
[1] 55450
ncol(vsdtrans_H3K27me3)
[1] 17
zscore_ATAC <- Xs_ATAC %>% as.data.frame() %>% tibble::rownames_to_column("ens_gene") %>% as_tibble
zscore_H3p3 <- Xs_H3p3 %>% as.data.frame() %>% tibble::rownames_to_column("ens_gene") %>% as_tibble
zscore_H3K27ac <- Xs_H3K27ac %>% as.data.frame() %>% tibble::rownames_to_column("ens_gene") %>% as_tibble
zscore_H3K4me3 <- Xs_H3K4me3 %>% as.data.frame() %>% tibble::rownames_to_column("ens_gene") %>% as_tibble
zscore_H3K27me3 <- Xs_H3K27me3 %>% as.data.frame() %>% tibble::rownames_to_column("ens_gene") %>% as_tibble
zscore <- zscore_ATAC %>% inner_join(zscore_H3p3) %>% inner_join(zscore_H3K27ac) %>% inner_join(zscore_H3K4me3) %>% inner_join(zscore_H3K27me3)
Joining, by = "ens_gene"
Joining, by = "ens_gene"
Joining, by = "ens_gene"
Joining, by = "ens_gene"
readr::write_csv(zscore, paste(csvfilepath,"_zscore_All.csv",sep="_"))
zscore_type <- zscore %>% inner_join(ensemble) %>% dplyr::select("ens_gene","ext_gene", "biotype","chr", all_of(def_list_select$sample))
Joining, by = "ens_gene"
normalized count, zscore の分布 (20191017修正)
norm_plotlist_all <- normalizedcount %>% gather("sample", "normalized",-(ens_gene)) %>% inner_join(def_list_select, by = "sample")
norm_plot_all <- norm_plotlist_all %>%
ggplot(aes(time_replicate,normalized,group=time_replicate,colour=type))+geom_violin(scale="count")+geom_boxplot(width=.1,)+facet_wrap(~seq*type,ncol=2)+theme_bw()+ theme(strip.text = element_text(size = 12),axis.text.x = element_text(hjust = 0.5,vjust = 0.5)) + ggtitle(paste(folder_name, "(all, normalized count)")) + scale_y_log10(limits = c(0.1,NA))
#+ggsci::scale_color_d3("category20")
print(norm_plot_all)
ggsave(plot=norm_plot_all,file="./NormCount.pdf", width = 20, height = 15, dpi = 360, limitsize = FALSE)
20191016修正, 20200623修正
chr_chart <- matome5 %>% dplyr::select(ens_gene,chr) %>% group_by(ens_gene,chr) %>% summarize() %>% ungroup()
`summarise()` regrouping output by 'ens_gene' (override with `.groups` argument)
#chr_chart <- bedfile %>% dplyr::select(ens_gene,chr) %>% group_by(ens_gene,chr) %>% summarize() %>% ungroup() #20191016修正
#> right_join(chr_chart, matome5, by="ens_gene")
# A tibble: 55,456 x 14
#bychr<- left_join(matome5, chr_chart, by="ens_gene") %>% dplyr::select(-(ens_gene)) %>%
# gather("sample","count",-chr) %>%
# group_by(chr,sample) %>% summarise(total=sum(count)) %>% ungroup
#mat <- left_join(matome5, chr_chart, by="ens_gene")
#bychr <- mat %>% dplyr::select(-(ens_gene)) %>%
# gather("sample","count",-chr) %>%
# group_by(chr,sample) %>% summarise(total=sum(count)) %>% f_sample() %>% ungroup
bychr <- matome5 %>% dplyr::select(chr,def_list_select$sample) %>%
gather("sample","count",-chr) %>%
group_by(chr,sample) %>% summarise(total=sum(count)) %>% f_sample() %>% ungroup
`summarise()` regrouping output by 'chr' (override with `.groups` argument)
ggplot(bychr,aes(reorder(sample,dplyr::desc(sample)),total/1e6,fill=chr)) +
theme_linedraw() + geom_bar(stat="identity") + coord_flip() +
xlab("sample") + ylab("million reads") + ggsci::scale_fill_igv() +
scale_x_discrete(limits = rev(levels(sample)))
NA
NA
NA
NA
mat_normcount <- left_join(normalizedcount, chr_chart, by="ens_gene")
bychr_mat_normcount <- mat_normcount %>% dplyr::select(-(ens_gene)) %>%
gather("sample","normcount",-chr) %>%
group_by(chr,sample) %>% summarise(normtotal=sum(normcount)) %>% f_sample() %>% ungroup
`summarise()` regrouping output by 'chr' (override with `.groups` argument)
ggplot(bychr_mat_normcount,aes(reorder(sample,dplyr::desc(sample)),normtotal/1e6,fill=chr)) +
theme_linedraw() + geom_bar(stat="identity") + coord_flip() +
xlab("sample") + ylab("normalized counts (million reads)") + ggsci::scale_fill_igv() +
scale_x_discrete(limits = rev(levels(sample)))
drop rows with all 0 -> +1/2 -> geom.scale -> log -> Pearson’s
#matf_Correlation <- mat %>% filter(chr!="chrM") %>% dplyr::select(-"chr") %>% filter_at(-(1),any_vars(. > 0))
matf_Correlation <- matome5 %>% dplyr::select(chr,def_list_select$sample) %>% filter(chr!="chrM") %>% dplyr::select(-"chr") %>% filter_at(-(1),any_vars(. > 0))
X_Correlation <- matf_Correlation %>% dplyr::select(-(1)) %>% as.matrix
rownames(X_Correlation) <- matf_Correlation$ens_gene
Unknown or uninitialised column: `ens_gene`.
lX_Correlation <- log(gscale(X_Correlation+0.5))
R <- cor(lX_Correlation); diag(R) <- NA
pheatmap::pheatmap(R,color=viridis::viridis(256))
#X <- matf %>% dplyr::select(-(1:4)) %>% as.matrix
#rownames(X) <- matf$ens_gene
#lX <- log(gscale(X+0.5))
#R <- cor(lX); diag(R) <- NA
#pheatmap::pheatmap(R,color=viridis::viridis(256))
###--- DESeq2によりnormalized count ---###
#model<- ~group
#dds_Correlation <- DESeq2::DESeqDataSetFromMatrix(X_Correlation[,def_list_select$sample],def_list_select,model) #必要なサンプルのみ
#dds_Correlation <- DESeq2::DESeqDataSetFromMatrix(X_Correlation[,def_list$sample],def_list,model)
#dds_Correlation <- DESeq2::estimateSizeFactors(dds_Correlation)
#norm_Correlation <- DESeq2::counts(dds_Correlation,normalized=TRUE) #normにnarmalized countが入る。
# set scale=TRUE if the patterns (not level) is the matter
p <- prcomp(t(lX_Correlation[rank(-apply(lX_Correlation,1,var)) <= ntop,]),scale=scalerows,center=TRUE)
screeplot(p,las=2,main="Importance")
print(summary(p)$imp[,seq(min(10,ncol(X_Correlation)))])
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Standard deviation 22.17018 1.317342 1.075728 0.783874 0.6637618 0.6246527 0.6047705 0.5830881 0.5684595 0.5323901
Proportion of Variance 0.98303 0.003470 0.002310 0.001230 0.0008800 0.0007800 0.0007300 0.0006800 0.0006500 0.0005700
Cumulative Proportion 0.98303 0.986500 0.988820 0.990050 0.9909300 0.9917100 0.9924400 0.9931200 0.9937700 0.9943300
label_Correlation <- def_list_select %>% filter(sample %in% colnames(X_Correlation))
df <- data.frame(p$x) %>% as_tibble(rownames="sample") %>%
inner_join(label_Correlation,.)
Joining, by = "sample"
print(df)
QC 終了
#------- setting -------#
fdr <- 0.1 # acceptable false discovery rate (固定)
lfcthreth <- log2(1) # threshold in abs(log2FC)
# controls should be placed in the right
plot_title1 <- "2gun"
contrast_H3p3 <- list(
#Intercept = list("Intercept"),
group_H3p3_UI_Doxplus_vs_minus = c("group", "H3p3_UI_DoxPlus","H3p3_UI_DoxMinus"),
group_H3p3_0h_Doxplus_vs_minus = c("group","H3p3_0h_DoxPlus", "H3p3_0h_DoxMinus"),
group_H3p3_24h_Doxplus_vs_minus = c("group","H3p3_24h_DoxPlus", "H3p3_24h_DoxMinus"),
group_H3p3_48h_Doxplus_vs_minus = c("group","H3p3_48h_DoxPlus", "H3p3_48h_DoxMinus")
)
contrast_H3K4me3 <- list(
#Intercept = list("Intercept"),
group_H3K4me3_UI_Doxplus_vs_minus = c("group", "H3K4me3_UI_DoxPlus","H3K4me3_UI_DoxMinus"),
group_H3K4me3_0h_Doxplus_vs_minus = c("group","H3K4me3_0h_DoxPlus", "H3K4me3_0h_DoxMinus"),
group_H3K4me3_24h_Doxplus_vs_minus = c("group","H3K4me3_24h_DoxPlus", "H3K4me3_24h_DoxMinus"),
group_H3K4me3_48h_Doxplus_vs_minus = c("group","H3K4me3_48h_DoxPlus", "H3K4me3_48h_DoxMinus")
)
contrast_H3K27ac <- list(
#Intercept = list("Intercept"),
group_H3K27ac_UI_Doxplus_vs_minus = c("group", "H3K27ac_UI_DoxPlus","H3K27ac_UI_DoxMinus"),
group_H3K27ac_0h_Doxplus_vs_minus = c("group","H3K27ac_0h_DoxPlus", "H3K27ac_0h_DoxMinus"),
group_H3K27ac_24h_Doxplus_vs_minus = c("group","H3K27ac_24h_DoxPlus", "H3K27ac_24h_DoxMinus"),
group_H3K27ac_48h_Doxplus_vs_minus = c("group","H3K27ac_48h_DoxPlus", "H3K27ac_48h_DoxMinus")
)
contrast_H3K27me3 <- list(
#Intercept = list("Intercept"),
group_H3K27me3_UI_Doxplus_vs_minus = c("group", "H3K27me3_UI_DoxPlus","H3K27me3_UI_DoxMinus"),
group_H3K27me3_0h_Doxplus_vs_minus = c("group","H3K27me3_0h_DoxPlus", "H3K27me3_0h_DoxMinus"),
group_H3K27me3_24h_Doxplus_vs_minus = c("group","H3K27me3_24h_DoxPlus", "H3K27me3_24h_DoxMinus"),
group_H3K27me3_48h_Doxplus_vs_minus = c("group","H3K27me3_48h_DoxPlus", "H3K27me3_48h_DoxMinus")
)
contrast_ATAC <- list(
#Intercept = list("Intercept"),
group_ATAC_UI_Doxplus_vs_minus = c("group", "ATAC_UI_DoxPlus","ATAC_UI_DoxMinus"),
group_ATAC_48h_Doxplus_vs_minus = c("group","ATAC_48h_DoxPlus", "ATAC_48h_DoxMinus")
)
# BRB
# group_UI_Doxplus_vs_minus = c("group", "BRB_UI_DoxPlus", "BRB_UI_DoxMinus"),
# group_0h_Doxplus_vs_minus = c("group", "BRB_0h_DoxPlus", "BRB_0h_DoxMinus"),
# group_24h_Doxplus_vs_minus = c("group", "BRB_24h_DoxPlus", "BRB_24h_DoxMinus"),
# group_48h_Doxplus_vs_minus = c("group", "BRB_48h_DoxPlus", "BRB_48h_DoxMinus")
#-----------------------#
#-------------- 自動 --------------------------------------------------#
#-- ファイル名 の設定 ---#
#folder_name_plot0 <- paste(".",folder_name, paste(folder_name,plot_title1,sep="_"),"",sep="/")
#folder_name_plot_path <- paste(folder_name_plot0,paste(folder_name,plot_title1,"",sep="_"),sep="")
log2 FC のみ計算 (DEGは特に出さない)
res_H3p3 <- mapply(function(x)
DESeq2::results(dds_H3p3,x,lfcThreshold=lfcthreth,alpha=fdr)
,contrast_H3p3)
res_H3K4me3 <- mapply(function(x)
DESeq2::results(dds_H3K4me3,x,lfcThreshold=lfcthreth,alpha=fdr)
,contrast_H3K4me3)
res_H3K27ac <- mapply(function(x)
DESeq2::results(dds_H3K27ac,x,lfcThreshold=lfcthreth,alpha=fdr)
,contrast_H3K27ac)
res_H3K27me3 <- mapply(function(x)
DESeq2::results(dds_H3K27me3,x,lfcThreshold=lfcthreth,alpha=fdr)
,contrast_H3K27me3)
res_ATAC <- mapply(function(x)
DESeq2::results(dds_ATAC,x,lfcThreshold=lfcthreth,alpha=fdr)
,contrast_ATAC)
print(fdr)
[1] 0.1
re_H3p3_all <- map(res_H3p3,as_tibble,rownames="ens_gene") %>%
tibble(aspect=factor(names(.),names(.)),data=.) %>%
mutate(data=map(data,annotate)) %>%
unnest(cols = "data")
re_H3K4me3_all <- map(res_H3K4me3,as_tibble,rownames="ens_gene") %>%
tibble(aspect=factor(names(.),names(.)),data=.) %>%
mutate(data=map(data,annotate)) %>%
unnest(cols = "data")
re_H3K27ac_all <- map(res_H3K27ac,as_tibble,rownames="ens_gene") %>%
tibble(aspect=factor(names(.),names(.)),data=.) %>%
mutate(data=map(data,annotate)) %>%
unnest(cols = "data")
re_H3K27me3_all <- map(res_H3K27me3,as_tibble,rownames="ens_gene") %>%
tibble(aspect=factor(names(.),names(.)),data=.) %>%
mutate(data=map(data,annotate)) %>%
unnest(cols = "data")
re_ATAC_all <- map(res_ATAC,as_tibble,rownames="ens_gene") %>%
tibble(aspect=factor(names(.),names(.)),data=.) %>%
mutate(data=map(data,annotate)) %>%
unnest(cols = "data")
filename <- paste("./2gun/",csvfilepath,"_resultsall_fdr0p1_H3p3.csv",sep="")
print(filename)
[1] "./2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3p3.csv"
readr::write_csv(re_H3p3_all,filename)
nrow(re_H3p3_all)
[1] 221800
filename <- paste("./2gun/",csvfilepath,"_resultsall_fdr0p1_H3K4me3.csv",sep="")
print(filename)
[1] "./2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3K4me3.csv"
readr::write_csv(re_H3K4me3_all,filename)
nrow(re_H3K4me3_all)
[1] 221800
filename <- paste("./2gun/",csvfilepath,"_resultsall_fdr0p1_H3K27ac.csv",sep="")
print(filename)
[1] "./2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3K27ac.csv"
readr::write_csv(re_H3K27ac_all,filename)
nrow(re_H3K27ac_all)
[1] 221800
filename <- paste("./2gun/",csvfilepath,"_resultsall_fdr0p1_H3K27me3.csv",sep="")
print(filename)
[1] "./2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_H3K27me3.csv"
readr::write_csv(re_H3K27me3_all,filename)
nrow(re_H3K27me3_all)
[1] 221800
filename <- paste("./2gun/",csvfilepath,"_resultsall_fdr0p1_ATAC.csv",sep="")
print(filename)
[1] "./2gun/TSS_pm5kb_20200624_resultsall_fdr0p1_ATAC.csv"
readr::write_csv(re_ATAC_all,filename)
nrow(re_ATAC_all)
[1] 110900
clustering H3.3
#20191205修正と作成
cluster_number <- 6
##--------- clustering -----------#
set.seed(3)
# H3.3で
zscore_H3p3_s <- zscore_type %>% dplyr::select("ens_gene",all_of(def_list_select_1$sample)) %>% filter(across(where(is_double), ~ (.x != 0)|(.x == 0)))
#zscore_H3p3_s <- zscore_H3p3 %>% filter(across(where(is_double), ~ (.x != 0)|(.x == 0)))
nrow(zscore_type)
[1] 55197
nrow(zscore_H3p3_s)
[1] 52326
Xs_H3p3 <- zscore_H3p3_s %>% dplyr::select(-ens_gene) %>% as.matrix()
rownames(Xs_H3p3) <- zscore_H3p3_s$ens_gene
km_allH3p3 <- kmeans(Xs_H3p3,cluster_number,nstart = 25,algorithm = "Lloyd")
10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした 10 回の反復を行いましたが収束しませんでした
kmc_allH3p3 <- km_allH3p3$centers %>% as_tibble(rownames="cluster") %>% gather(sample,val,-cluster) %>% inner_join(def_list_select)
Joining, by = "sample"
kmc_allH3p3_group <- kmc_allH3p3
#kmc_LRT_group <- kmc_LRT %>% mutate(growth=factor(growth, c("UI","Diff0h","Diff24h","Diff48h"))) %>% mutate(type=factor(type, c("Doxplus","Doxminus")))
#kmc_LRT_group <- kmc_LRT_group %>% mutate(time=case_when(growth=="UI" ~"UI",growth=="Diff0h"~"0h",growth=="Diff24h"~"24h",growth=="Diff48h"~"48h",TRUE~"error"))
#kmc_LRT_group <- kmc_LRT_group %>% mutate(time=factor(time, c("UI","0h","24h","48h")))
#gggglabel <- paste("k-means: Total",nrow(Xs_H3p3),"[1]",km_allH3p3$size[1],"[2]",km_allH3p3$size[2],"[3]",km_allH3p3$size[3],"[4]",km_allH3p3$size[4],"[5]",km_allH3p3$size[5],"[6]",km_allH3p3$size[6],sep=" ")
gggglabel <- paste("Original",nrow(zscore_type),"H3.3 k-means: Total",nrow(Xs_H3p3),"[1]",km_allH3p3$size[1],"[2]",km_allH3p3$size[2],"[3]",km_allH3p3$size[3],"[4]",km_allH3p3$size[4],"[5]",km_allH3p3$size[5],"[6]",km_allH3p3$size[6],sep=" ")
#------- size -------#
print(km_allH3p3$size)
[1] 8017 8143 10775 8039 7716 9636
#rrres_allH3p3 <- km_allH3p3$cluster %>% tibble(ens_gene=names(.),cluster=.) %>% left_join(zscore_type_clus3,.) %>% arrange(cluster) %>% dplyr::select(ens_gene,ext_gene,biotype,chr,cluster,all_of(def_list_select$sample))
rrres_allH3p3 <- km_allH3p3$cluster %>% tibble(ens_gene=names(.),cluster=.) %>% left_join(.,zscore_type) %>% arrange(cluster) %>% dplyr::select(ens_gene,ext_gene,biotype,chr,cluster,all_of(def_list_select$sample))
Joining, by = "ens_gene"
#rrres_allH3p3 <- km_allH3p3$cluster %>% tibble(ens_gene=names(.),cluster=.) %>% left_join(zscore_type,.) %>% arrange(cluster) %>% dplyr::select(ens_gene,ext_gene,biotype,chr,cluster,all_of(def_list_select$sample))
#rrres_LRT <- km_LRT$cluster %>% tibble(ens_gene=names(.),cluster=.) %>% right_join(e2g,.) %>% arrange(cluster)
file_path <- paste("./H3p3allcluster/", csvfilepath, "_kmeans_cluster.csv",sep="")
readr::write_csv(rrres_allH3p3,file_path)
##------- PCA -------#
pcacluster_save <- prcomp(Xs_H3p3)$x %>% as_tibble %>% dplyr::select(PC1,PC2) %>% mutate(cluster=km_allH3p3$cluster) %>% ggplot(aes(PC1,PC2,colour=factor(cluster)))+geom_point(size=1.5,alpha=0.6)+coord_fixed()+theme_linedraw()+ggsci::scale_color_d3("category20")
file_path <- paste("./H3p3allcluster/", csvfilepath, "_kmeans__pcacluster_PC1PC2.pdf",sep="")
ggsave(plot=pcacluster_save,file=file_path, width = 10, height = 6, dpi = 120)
print(pcacluster_save)
pcacluster_save <- prcomp(Xs_H3p3)$x %>% as_tibble %>% dplyr::select(PC1,PC3) %>% mutate(cluster=km_allH3p3$cluster) %>% ggplot(aes(PC1,PC3,colour=factor(cluster)))+geom_point(size=1.5,alpha=0.6)+coord_fixed()+theme_linedraw()+ggsci::scale_color_d3("category20")
file_path <- paste("./H3p3allcluster/", csvfilepath, "_kmeans__pcacluster_PC1PC3.pdf",sep="")
ggsave(plot=pcacluster_save,file=file_path, width = 10, height = 6, dpi = 120)
print(pcacluster_save)
#================================================#
# mouseCTX 0438を参考に。
#------------------#
f_cluster <- function(x) x %>% group_by(group, type, time, cluster, seq) %>% summarise(avg=mean(val),se=sd(val)/sqrt(length(val))) %>% ungroup()
print(kmc_allH3p3_group %>% group_by(group, type, time) %>% summarise())
`summarise()` regrouping output by 'group', 'type' (override with `.groups` argument)
f_clusterp <- function(x) x %>% group_by(group, type, time, cluster, sep) %>% summarise(avg=mean(val),se=sd(val)/sqrt(length(val))) %>% ungroup()
print(kmc_allH3p3_group %>% group_by(group, type, time) %>% summarise()) #作図用
`summarise()` regrouping output by 'group', 'type' (override with `.groups` argument)
#-------#
cluster_save <- kmc_allH3p3_group %>%
ggplot(aes(time,val,group=type,colour=type))+ geom_abline(intercept=0,slope=0,linetype="dashed",colour="gray") +geom_line(aes(x=time,y=avg,colour=type),data=f_cluster)+geom_point()+facet_wrap(~cluster*seq,ncol=3)+ggsci::scale_color_npg()+theme_bw()+ theme(strip.text = element_text(size = 12),axis.text.x = element_text(vjust = 0.5), strip.background = element_blank(), plot.title=element_text(size=5)) + ggtitle(gggglabel)+ggsci::scale_color_npg() + ylab("z score")
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
file_path <- paste("./H3p3allcluster/", csvfilepath, "_cluster_type.pdf",sep="")
ggsave(plot=cluster_save,file=file_path, width = 6, height = 6, dpi = 120)
#ggsave(plot=cluster_save,file=file_path, width = 6, height = 6, dpi = 120)
print(cluster_save)
#================================================#
z_H3p3clus1 <- rrres_allH3p3 %>% filter(cluster=="1") %>% dplyr::rename(H3p3clus=cluster)
z_H3p3clus2 <- rrres_allH3p3 %>% filter(cluster=="2") %>% dplyr::rename(H3p3clus=cluster)
z_H3p3clus3 <- rrres_allH3p3 %>% filter(cluster=="3") %>% dplyr::rename(H3p3clus=cluster)
z_H3p3clus4 <- rrres_allH3p3 %>% filter(cluster=="4") %>% dplyr::rename(H3p3clus=cluster)
z_H3p3clus5 <- rrres_allH3p3 %>% filter(cluster=="5") %>% dplyr::rename(H3p3clus=cluster)
z_H3p3clus6 <- rrres_allH3p3 %>% filter(cluster=="6") %>% dplyr::rename(H3p3clus=cluster)
nrow(rrres_allH3p3)
[1] 52326
nrow(z_H3p3clus1)
[1] 8017
nrow(z_H3p3clus2)
[1] 8143
nrow(z_H3p3clus3)
[1] 10775
nrow(z_H3p3clus4)
[1] 8039
nrow(z_H3p3clus5)
[1] 7716
nrow(z_H3p3clus6)
[1] 9636
nrow(rrres_allH3p3 %>% filter(is.na(cluster)))
[1] 0
rrres_allH3p3 %>% filter(ext_gene %in% c("Myh3","Ckm","Acta1","Tnnt2","Actb","Csrp3","Tpm2","Nsdhl","Myog"))
cluster_num <- 4
##---- BRBのクラスタリング(LRT) の結果 -------#
#filepath_BRBcluster <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/LRT/DEG_fdr0p1__BRB0432lane2noumi_H3mm18_Dox_kmeans4__cluster_result.csv" #BRB等のDEGのリスト
#filepath_BRBallgene <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/H3mm18KO_3T3_Dox_normCount_genename.csv" #BRB等のDEGのリスト
#filepath_BRBzscore <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/LRT/clustering_XsLRTall__BRB0432lane2noumi_H3mm18_Dox.csv"
#filepath_BRBdef <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/deftable_BRB_noumi_new_190520_fin191205ver.txt"
#filepath_BRBdef <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/deftable_BRB_noumi_new_190520_fin191205ver_20200625.txt"
#filepath_BRBChILATAC_def <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/TSS_count/ChILAll_TSS_pm5kb_withATAC/deftable_multicov_ChIL01100111_20200501_3T3_EGFP18_UI_DoxMinus_H3p3K27acK4Kme327me3_withATAC_withBRB.txt"
#filepath_BRBensemble <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/ensemble_list_asia.csv" #BRBの時のgene名リスト
#filepath_BRB_FC_all <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/LRT/all__BRB0432lane2noumi_H3mm18_Dox.csv" #BRB等のDEGのリスト
#filepath_BRB_2gun <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/2gun/BRB0432lane2noumi_H3mm18_Dox_results_fdr0p1__final191205.csv" #BRB等のDEGのリスト
filepath_BRBdef <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/deftable_BRB_noumi_new_190520_Last20200811ver.txt"
filepath_BRBcluster <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/LRT/DEG_fdr0p1__BRB0432lane2noumi_H3mm18_Dox_kmeans4__cluster_result.csv" #BRB等のDEGのリスト
filepath_BRBallgene <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/H3mm18KO_3T3_Dox_normCount_genename.csv" #BRB等のDEGのリスト
filepath_BRBzscore <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/H3mm18KO_3T3_Dox__zscore_type_all.csv"
##20200811 log2FC
filepath_BRB_FC_deseq <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/2gun/BRB0432lane2noumi_H3mm18_Dox_resultsall_fdr0p1__final191205_last200811.csv" #BRB等のDEGのリスト
filepath_BRB_normcount <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Last_Rserver_200811/BRB0432lane2noumi_H3mm18_Dox__normCount__final191205_last200811.csv"
#filepath_BRB_normcount <- "/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/BRB0432lane2noumi_H3mm18_Dox__normCount__final191205.csv"
#Final_Last_Rserver_200811
BRBのデータ (DEG,log2FCを抜き出す)
#--- DEGのリストの呼び出し --------#
cluster_BRBlist <- readr::read_csv(filepath_BRBcluster) %>% mutate(cluster=factor(cluster,c(1:cluster_num))) #BRB等のDEGのリスト
Parsed with column specification:
cols(
ens_gene = col_character(),
ext_gene = col_character(),
biotype = col_character(),
chr = col_character(),
cluster = col_double(),
baseMean = col_double(),
log2FoldChange = col_double(),
lfcSE = col_double(),
stat = col_double(),
pvalue = col_double(),
padj = col_double()
)
#--- log2 FCのリストの呼び出し --------#
re_BRB_all <- readr::read_csv(filepath_BRB_FC_deseq)
Parsed with column specification:
cols(
aspect = col_character(),
ens_gene = col_character(),
ext_gene = col_character(),
biotype = col_character(),
chr = col_character(),
baseMean = col_double(),
log2FoldChange = col_double(),
lfcSE = col_double(),
stat = col_double(),
pvalue = col_double(),
padj = col_double()
)
#--- BRBの全geneリストの呼び出し --------#
all_BRBlist <- readr::read_csv(filepath_BRBallgene) #normcount
Parsed with column specification:
cols(
.default = col_double(),
ens_gene = col_character(),
ext_gene = col_character(),
biotype = col_character(),
chr = col_character()
)
See spec(...) for full column specifications.
nrow(all_BRBlist)
[1] 21707
#----#
zscore_BRB <- readr::read_csv(filepath_BRBzscore)
Parsed with column specification:
cols(
.default = col_double(),
ens_gene = col_character(),
ext_gene = col_character(),
biotype = col_character(),
chr = col_character()
)
See spec(...) for full column specifications.
nrow(zscore_BRB)
[1] 21707
#----#
def_BRB <- readr::read_tsv(filepath_BRBdef) %>% mutate(seq=factor(seq, c("BRB"))) %>%
mutate(time1 = time, rep1=rep) %>% unite(time1,rep1,col="time_replicate") %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% mutate(type=factor(type,c("DoxPlus","DoxMinus"))) %>% mutate(rep=factor(rep, c("1", "2", "3", "4")))%>% mutate(time_replicate=factor(time_replicate,c("UI_1", "UI_2", "UI_3", "UI_4", "0h_1","0h_2","24h_1","24h_2","48h_1","48h_2","48h_3","48h_4")))
Parsed with column specification:
cols(
file = col_character(),
sample0 = col_character(),
barcode = col_character(),
growth = col_character(),
sample = col_character(),
group = col_character(),
time = col_character(),
type = col_character(),
seq = col_character(),
rep = col_double()
)
#----#
#FC_rank_all_BRBlist <- re_BRB_all %>% arrange(desc(abs(log2FoldChange))) %>% mutate(Rank=row_number()) #normcount
#nrow(FC_rank_all_BRBlist)
#print(FC_rank_all_BRBlist)
#%>% mutate(cluster=factor(cluster,c(1:cluster_num))) #BRB等のDEGのリスト
###--- DESeq2により計算したnormalized countのうち、BRB等のDEGのリストにあったものを書き出し ---###
#norm_table_select <- normalizedcount %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% full_join(cluster_BRBlist, ., by="ens_gene") #DEGリストのnormalized count
#readr::write_csv(norm_table_select, paste(".",folder_name, paste(RNAseq_cluster,"__normalizedcount__",folder_name,".csv",sep=""),sep="/")) #normalized count csvにgene名を付けた。
#------------------------------------------------------------------------------------------------#
#----#
norm_BRB_def_original <- readr::read_csv(filepath_BRB_normcount)
Parsed with column specification:
cols(
ens_gene = col_character(),
sample = col_character(),
normalized = col_double(),
file = col_character(),
sample0 = col_character(),
barcode = col_character(),
growth = col_character(),
group = col_character(),
time = col_character(),
type = col_character(),
seq = col_character(),
rep = col_double(),
count = col_double(),
ext_gene = col_character(),
biotype = col_character(),
chr = col_character()
)
#norm_BRB <- norm_BRB_def_original %>% dplyr::rename(Type=type,rep=replicate,norm=normalized) %>% mutate(type=case_when(Type=="Doxminus"~"DoxMinus",Type=="Doxplus"~"DoxPlus")) %>% mutate(seq="BRB") %>%
#mutate(time1 = time, rep1=rep) %>% unite(time1,rep1,col="time_replicate") %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% mutate(type=factor(type,c("DoxPlus","DoxMinus"))) %>% mutate(rep=factor(rep, c("1", "2", "3", "4")))%>% mutate(time_replicate=factor(time_replicate,c("UI_1", "UI_2", "UI_3", "UI_4", "0h_1","0h_2","24h_1","24h_2","48h_1","48h_2","48h_3","48h_4")))
norm_BRB <- norm_BRB_def_original %>% dplyr::rename(norm=normalized) %>% mutate(seq="BRB") %>%
mutate(time1 = time, rep1=rep) %>% unite(time1,rep1,col="time_replicate") %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% mutate(type=factor(type,c("DoxPlus","DoxMinus"))) %>% mutate(rep=factor(rep, c("1", "2", "3", "4")))%>% mutate(time_replicate=factor(time_replicate,c("UI_1", "UI_2", "UI_3", "UI_4", "0h_1","0h_2","24h_1","24h_2","48h_1","48h_2","48h_3","48h_4")))
# %>% mutate(group=case_when(Group=="Doxminus_Diff0h"~"BRB_0h_DoxMinus",Group=="Doxplus_Diff0h"~"BRB_0h_DoxPlus"))
print(norm_BRB)
#----#
#BRB_2gun <- readr::read_csv(filepath_BRB_2gun)
#"/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/LRT/clustering_XsLRTall__BRB0432lane2noumi_H3mm18_Dox.csv"
#--- DEGのリストに位置情報を加える
#bedfile_cluster <- bedfile %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% dplyr::select(ens_gene,TSS_region,gene_region,Strand) %>% full_join(cluster_BRBlist, ., by="ens_gene") #位置情報あり
# dplyr::select(ens_gene,TSS_region,gene_region,Strand)
#NormCountBRBmat <- readr::read_csv("/home/guestA/o70578a/akuwakado/kuwakado/BRBSeq/H3mm18_Dox_0432lane2/Final_Rserver_191203/H3mm18KO_3T3_Dox_normCount_genename.csv")
list1 <- rrres_allH3p3 %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% group_by(cluster) %>% summarise(BRBcount=n(),genelist=paste(ext_gene,collapse=", "),IDlist=paste(ens_gene,collapse=", "))
`summarise()` ungrouping output (override with `.groups` argument)
#list2 <- rrres_H3p3clus2 %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% group_by(cluster) %>% summarise(BRBcount=n(),genelist=paste(ext_gene,collapse=", "),IDlist=paste(ens_gene,collapse=", "))
readr::write_csv(list1,"Compare_ChILATAC_cluster_vs_BRBDEGs.csv")
#readr::write_csv(list2,"Compare_ChILATAC_cluster2_vs_BRBDEGs.csv")
print(list1)
cluster1 <- cluster_BRBlist %>% filter(cluster=="1")
cluster2 <- cluster_BRBlist %>% filter(cluster=="2")
cluster3 <- cluster_BRBlist %>% filter(cluster=="3")
cluster4 <- cluster_BRBlist %>% filter(cluster=="4")
listclus <- rrres_allH3p3 %>% dplyr::select(ens_gene,ext_gene,cluster) %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% left_join(dplyr::select(cluster_BRBlist,ens_gene,cluster) %>% dplyr::rename(BRBclus=cluster)) %>% group_by(cluster,BRBclus) %>% summarise(count=n(),genelist=paste(ext_gene,collapse=", "),IDlist=paste(ens_gene,collapse=", "))
Joining, by = "ens_gene"
`summarise()` regrouping output by 'cluster' (override with `.groups` argument)
listclus2 <- listclus %>% dplyr::select(cluster,BRBclus,count) %>% mutate(BRBclus=paste("BRBcluster",BRBclus,sep="")) %>% spread(key=BRBclus,value=count, fill = 0)
readr::write_csv(listclus,"./Compare/Compare_ChILATAC_cluster_vs_BRBDEGs_summary.csv")
readr::write_csv(listclus2,"./Compare/Compare_ChILATAC_cluster_vs_BRBDEGs_summarycount.csv")
print(listclus)
print(listclus2)
NA
NA
#---#
H3p3BRBlistclus <- listclus2 %>% mutate(H3p3cluster=paste(cluster,sep = ""))
#H3p3BRBlistclus <- listclus2 %>% mutate(H3p3cluster=paste("H3p3cluster",cluster,sep = ""))
H3p3BRBlistclus_mat <- H3p3BRBlistclus %>% ungroup() %>% dplyr::select(-cluster,-H3p3cluster) %>% as.matrix()
rownames(H3p3BRBlistclus_mat) <- H3p3BRBlistclus$H3p3cluster
H3p3BRBlistclus_count <- H3p3BRBlistclus %>% ungroup() %>% dplyr::select(-cluster) %>% gather(key=sample,value=count,-H3p3cluster)
#---#
resultchisq <- chisq.test(H3p3BRBlistclus_mat)
カイ自乗近似は不正確かもしれません
resultchisq
Pearson's Chi-squared test
data: H3p3BRBlistclus_mat
X-squared = 37.639, df = 15, p-value = 0.00102
#resultchisq$residuals
#resultchisq$expected
resultchisq$expected %>% sum()
[1] 226
#residuals <- resultchisq$residuals %>% as.data.frame(.) %>% tibble::rownames_to_column("H3p3cluster") %>% gather(key=sample,value=value,-(H3p3cluster)) %>% mutate(H3p3cluster=factor(H3p3cluster,c("H3p3cluster6","H3p3cluster5","H3p3cluster4","H3p3cluster3","H3p3cluster2","H3p3cluster1")))
#H3p3BRBlistclus_count_residuals <- residuals %>% left_join(H3p3BRBlistclus_count) %>% mutate(H3p3cluster=factor(H3p3cluster,c("H3p3cluster6","H3p3cluster5","H3p3cluster4","H3p3cluster3","H3p3cluster2","H3p3cluster1")))
residuals <- resultchisq$residuals %>% as.data.frame(.) %>% tibble::rownames_to_column("H3p3cluster") %>% gather(key=sample,value=value,-(H3p3cluster)) %>% mutate(H3p3cluster=factor(H3p3cluster,c("1","2","3","4","5","6"))) %>% mutate(BRBcluster=gsub("BRBcluster","",sample)) %>% mutate(BRBcluster=factor(BRBcluster,c("4","3","2","1")))
H3p3BRBlistclus_count_residuals <- residuals %>% left_join(H3p3BRBlistclus_count)
Joining, by = c("H3p3cluster", "sample")
#%>% mutate(H3p3cluster=factor(H3p3cluster,c("1","2","3","4","5","6")))
paste(resultchisq$method,"Residual",sep=": ")
[1] "Pearson's Chi-squared test: Residual"
chisq_plot <- H3p3BRBlistclus_count_residuals %>% ggplot(aes(x=H3p3cluster, y=BRBcluster, fill=value, label=count)) + geom_tile() + geom_text(aes(x=H3p3cluster, y=BRBcluster, label=as.character(count))) + scale_fill_gradient2(low="blue", high="red", na.value="black", name="") + theme(axis.text.x = element_text(angle = 90),title = element_text(size=2),legend.position = "top") + theme_minimal() + ylab("BRB DEG cluster") + xlab("H3.3 cluster")
chisq_plot
ggsave(file="./Compare/H3p3BRBlistclus.pdf", plot = chisq_plot, dpi = 100, width = 5, height = 3,limitsize = FALSE)
chisq_plot2 <- H3p3BRBlistclus_count_residuals %>% ggplot(aes(x=H3p3cluster, y=BRBcluster, fill=value, label=count)) + geom_tile() + scale_fill_gradient2(low="blue", high="red", na.value="black", name="") + theme(axis.text.x = element_text(angle = 90),title = element_text(size=2),legend.position = "top") + theme_minimal() + ylab("BRB DEG cluster") + xlab("H3.3 cluster")
chisq_plot2
ggsave(file="./Compare/H3p3BRBlistclus_notitle.pdf", plot = chisq_plot2, dpi = 100, width = 5, height = 3,limitsize = FALSE)
H3p3BRBlistclus_count_residuals %>% readr::write_csv("./Compare/H3p3BRBlistclus.csv")
##——— リストを保存 ————-# #– 確認 –#
rrres_allH3p3 %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% left_join(cluster_BRBlist %>% dplyr::rename(BRBclus=cluster)) %>% dplyr::select(cluster)%>% group_by(cluster) %>% summarise(count=n())
rrres_allH3p3 %>% filter(ens_gene %in% cluster_BRBlist$ens_gene) %>% left_join(cluster_BRBlist %>% dplyr::rename(BRBclus=cluster)) %>% dplyr::select(BRBclus,cluster) %>% group_by(BRBclus,cluster) %>% summarise(count=n())
##————————————#
normalized count (BRB DEG) の deftable等 (発現が低いものはCut)
まず、BRBで遺伝子発現が十分大きいものに絞り込む(Cut off リストの作成)
Set_cutoff <- 10.0
## 各時刻の平均を計算し、normalized count > 10 を超えるものを抽出する。
#----- SKMとCTXのみ取り出す ---# 20191205
#norm_BRB_all <- norm_BRB %>% gather("sample", "normalized",-(ens_gene)) %>% inner_join(def, by = "sample")
#norm_BRB_all <- norm_BRB_all %>% filter(intact_CTX=="CTX"|intact_CTX=="SKM") %>% mutate(WT_KO=factor(WT_KO, c("H3mm18KO","WT"))) %>% mutate(Day=factor(Day, c("Day0","Day5","Day14"))) %>% mutate(intact_CTX=factor(intact_CTX, c("CTX","SKM")))
#notm_plotlist_cutoff <- norm_plotlist_all %>% annotate() %>% group_by(ens_gene, ext_gene, Day, intact_CTX) %>% summarize(groupMean=mean(normalized)) %>% ungroup() %>% dplyr::select(ens_gene, ext_gene) %>% unique()
norm_BRB_beforecutoff <- norm_BRB %>% group_by(ens_gene, ext_gene, seq, time) %>% summarize(groupMean=mean(norm))
`summarise()` regrouping output by 'ens_gene', 'ext_gene', 'seq' (override with `.groups` argument)
nrow(norm_BRB_beforecutoff)
[1] 86828
nrow(norm_BRB_beforecutoff %>% ungroup() %>% dplyr::select(ens_gene, ext_gene) %>% unique()) #この値をMAplotのx軸に使用
[1] 21707
print("--- cut off ---")
[1] "--- cut off ---"
norm_BRB_cutoff <- norm_BRB_beforecutoff %>% filter(groupMean > Set_cutoff) %>% ungroup()
nrow(norm_BRB_cutoff)
[1] 28442
norm_BRB_cutoff_list <-norm_BRB_cutoff %>% dplyr::select(ens_gene, ext_gene) %>% unique()
nrow(norm_BRB_cutoff_list)
[1] 8159
norm_BRB_beforecutoff %>% readr::write_csv("./log2FC/tables/Norm_BRB_groupMean.csv")
norm_BRB_cutoff %>% readr::write_csv("log2FC/tables/Norm_BRB_groupMean_cutoff10.csv")
norm_BRB_cutoff_list %>% readr::write_csv("log2FC/tables/Norm_BRB_groupMean_cutoff10_genelist.csv")
re_H3p3_FC_cutoff <- re_H3p3_all %>% mutate(seq="H3p3", time=gsub("group_H3p3_","",aspect)) %>% mutate(time=gsub("_Doxplus_vs_minus","",time)) %>% dplyr::select(ens_gene,ext_gene,biotype,chr, log2FoldChange,seq,time) %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene)
re_H3K4me3_FC_cutoff <- re_H3K4me3_all %>% mutate(seq="H3K4me3", time=gsub("group_H3K4me3_","",aspect)) %>% mutate(time=gsub("_Doxplus_vs_minus","",time)) %>% dplyr::select(ens_gene,ext_gene,biotype,chr, log2FoldChange,seq,time) %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene)
re_H3K27ac_FC_cutoff <- re_H3K27ac_all %>% mutate(seq="H3K27ac", time=gsub("group_H3K27ac_","",aspect)) %>% mutate(time=gsub("_Doxplus_vs_minus","",time)) %>% dplyr::select(ens_gene,ext_gene,biotype,chr, log2FoldChange,seq,time) %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene)
re_H3K27me3_FC_cutoff <- re_H3K27me3_all %>% mutate(seq="H3K27me3", time=gsub("group_H3K27me3_","",aspect)) %>% mutate(time=gsub("_Doxplus_vs_minus","",time)) %>% dplyr::select(ens_gene,ext_gene,biotype,chr, log2FoldChange,seq,time) %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene)
re_ATAC_FC_cutoff <- re_ATAC_all %>% mutate(seq="ATAC", time=gsub("group_ATAC_","",aspect)) %>% mutate(time=gsub("_Doxplus_vs_minus","",time)) %>% dplyr::select(ens_gene,ext_gene,biotype,chr, log2FoldChange,seq,time) %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene)
re_BRB_FC_cutoff <- re_BRB_all %>% mutate(seq="BRB", time=gsub("group_","",aspect)) %>% mutate(time=gsub("_Doxplus_vs_minus","",time)) %>% dplyr::select(ens_gene,ext_gene,biotype,chr, log2FoldChange,seq,time) %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene)
filename <- "./log2FC/tables/log2FC_H3p3_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/log2FC_H3p3_cutoff10.csv"
readr::write_csv(re_H3p3_FC_cutoff,filename)
nrow(re_H3p3_FC_cutoff)
[1] 32636
filename <- "./log2FC/tables/log2FC_H3K4me3_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/log2FC_H3K4me3_cutoff10.csv"
readr::write_csv(re_H3K4me3_FC_cutoff,filename)
nrow(re_H3K4me3_FC_cutoff)
[1] 32636
filename <- "./log2FC/tables/log2FC_H3K27ac_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/log2FC_H3K27ac_cutoff10.csv"
readr::write_csv(re_H3K27ac_FC_cutoff,filename)
nrow(re_H3K27ac_FC_cutoff)
[1] 32636
filename <- "./log2FC/tables/log2FC_H3K27me3_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/log2FC_H3K27me3_cutoff10.csv"
readr::write_csv(re_H3K27me3_FC_cutoff,filename)
nrow(re_H3K27me3_FC_cutoff)
[1] 32636
filename <- "./log2FC/tables/lo2gFC_ATAC_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/lo2gFC_ATAC_cutoff10.csv"
readr::write_csv(re_ATAC_FC_cutoff,filename)
nrow(re_ATAC_FC_cutoff)
[1] 16318
filename <- "./log2FC/tables/log2FC_BRB_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/log2FC_BRB_cutoff10.csv"
readr::write_csv(re_BRB_FC_cutoff,filename)
nrow(re_BRB_FC_cutoff)
[1] 32636
## 全て結合
re_all_FC_cutoff <- bind_rows(re_H3p3_FC_cutoff, re_H3K4me3_FC_cutoff) %>% bind_rows(re_H3K27ac_FC_cutoff) %>% bind_rows(re_H3K27me3_FC_cutoff) %>% bind_rows(re_ATAC_FC_cutoff) %>% bind_rows(re_BRB_FC_cutoff) %>% mutate(seq=factor(seq, c("H3p3", "H3K4me3","H3K27ac","H3K27me3","ATAC","BRB"))) %>% mutate(time=factor(time, c("UI", "0h","24h","48h")))
re_all_FC_cutoff
filename <- "./log2FC/tables/log2FC_ChILATACBRB_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/log2FC_ChILATACBRB_cutoff10.csv"
readr::write_csv(re_all_FC_cutoff,filename)
nrow(re_all_FC_cutoff)
[1] 179498
spread_all_FC_cutoff <- re_all_FC_cutoff %>% group_by(ens_gene, ext_gene, biotype, chr, time) %>% spread(key=seq,value=log2FoldChange)
nrow(spread_all_FC_cutoff)
[1] 32636
spread_all_FC_cutoff <- spread_all_FC_cutoff %>% left_join(dplyr::select(norm_BRB_cutoff, ens_gene, ext_gene, time, groupMean) %>% dplyr::rename(BRBgroupMean=groupMean))
Joining, by = c("ens_gene", "ext_gene", "time")
nrow(spread_all_FC_cutoff)
[1] 32636
filename <- "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10.csv"
print(filename)
[1] "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10.csv"
readr::write_csv(spread_all_FC_cutoff,filename)
spread_all_FC_cutoff_clus <- spread_all_FC_cutoff %>% left_join(dplyr::select(rrres_allH3p3,ens_gene,cluster)) %>% dplyr::rename(H3p3cluster=cluster) %>% left_join(dplyr::select(cluster_BRBlist,ens_gene,cluster)) %>% dplyr::rename(BRBDEGcluster=cluster)
Joining, by = "ens_gene"
Joining, by = "ens_gene"
filename <- "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10__withCluster.csv"
print(filename)
[1] "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10__withCluster.csv"
readr::write_csv(spread_all_FC_cutoff_clus,filename)
nrow(spread_all_FC_cutoff_clus)
[1] 32636
spread_all_FC_cutoff_clus %>% ungroup() %>% dplyr::select(ens_gene,H3p3cluster) %>% unique() %>% group_by(H3p3cluster) %>% summarize(H3p3_cutoff_count=n())
`summarise()` ungrouping output (override with `.groups` argument)
rrres_allH3p3 %>% group_by(cluster) %>% summarize(H3p3_geneAllcount=n())
`summarise()` ungrouping output (override with `.groups` argument)
f_gene_H3p3clus3 <- function(x) x %>% filter(H3p3cluster=="3")
f_gene_BRBclus3 <- function(x) x %>% filter(BRBDEGcluster=="3")
list_gene_qpcr <- c("Acta1","Myh3","Ttn","Myog")
spread_all_FC_cutoff_H3p3clus3 <- spread_all_FC_cutoff_clus %>% ungroup() %>% f_gene_H3p3clus3
nrow(spread_all_FC_cutoff_H3p3clus3)
[1] 13452
filename <- "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10__withCluster__H3p3clus3.csv"
print(filename)
[1] "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10__withCluster__H3p3clus3.csv"
readr::write_csv(spread_all_FC_cutoff_H3p3clus3,filename)
spread_all_FC_cutoff_H3p3clus3
spread_all_FC_cutoff_H3p3clus3 %>% f_gene_H3p3clus3 %>% f_gene_BRBclus3
spread_all_FC_cutoff_H3p3clus3 %>% f_gene_H3p3clus3 %>% f_gene_BRBclus3 %>% filter(ext_gene %in% list_gene_qpcr)
H3p3clus3cutoff <- spread_all_FC_cutoff_H3p3clus3 %>% ungroup() %>% dplyr::select(ens_gene) %>% unique() %>% nrow()
H3p3clus3cutoff_brbclus3 <- spread_all_FC_cutoff_H3p3clus3 %>% f_gene_BRBclus3 %>% ungroup() %>% dplyr::select(ens_gene) %>% unique() %>% nrow()
nrow(z_H3p3clus3)
[1] 10775
H3p3clus3cutoff
[1] 3363
H3p3clus3cutoff_brbclus3
[1] 55
plot_all_FC_cutoff_H3p3clus3 <- spread_all_FC_cutoff_H3p3clus3 %>% dplyr::mutate(label_text = dplyr::case_when(ext_gene %in% list_gene_qpcr ~ ext_gene, TRUE ~ ""),shape = dplyr::case_when(ext_gene %in% list_gene_qpcr ~ "TRUE", TRUE ~ "FALSE"))
nrow(plot_all_FC_cutoff_H3p3clus3)
[1] 13452
plot_all_FC_cutoff_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(!is.na(BRBgroupMean)) # groupMean > 10のみ残す
nrow(plot_all_FC_cutoff_H3p3clus3)
[1] 11802
filename <- "./log2FC/tables/Plot_log2FC_ChILATACBRB_cutoff10__withCluster__H3p3clus3.csv"
print(filename)
[1] "./log2FC/tables/Plot_log2FC_ChILATACBRB_cutoff10__withCluster__H3p3clus3.csv"
readr::write_csv(plot_all_FC_cutoff_H3p3clus3,filename)
plot_all_FC_cutoff_H3p3clus3 %>% group_by(time) %>% summarise(count=n()) #図中の数
`summarise()` ungrouping output (override with `.groups` argument)
plot_all_FC_cutoff_H3p3clus3 %>% f_gene_BRBclus3 %>% group_by(time) %>% summarise(count=n()) #図中の数(BRB DEG cluster3のみ)
`summarise()` ungrouping output (override with `.groups` argument)
正規分布ならピアソンだが、今回正規分布ではないのでスピアマンの順位相関係数 を使う
Time_list <- c("UI","0h","24h","48h")
Count_FC_cutoff_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% ungroup() %>% group_by(time) %>% summarise(Plot_genes=n()) #図中の数
`summarise()` ungrouping output (override with `.groups` argument)
#######
print("~~ H3p3_BRB ~~")
[1] "~~ H3p3_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clusterAll: H3p3_BRB --"))
corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3$H3p3, corr_H3p3clus3$BRB, method="spearman")
#print(tttttt)
if (i == 1) {
cortest_result_s <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("H3p3", "BRB",sep="_"))
}
else {
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("H3p3", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
}
[1] "----- UI --- H3p3clusterAll: H3p3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 0h --- H3p3clusterAll: H3p3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 24h --- H3p3clusterAll: H3p3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 48h --- H3p3clusterAll: H3p3_BRB --"
タイのため正確な p 値を計算することができません
print("~~ H3K4me3_BRB ~~")
[1] "~~ H3K4me3_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clusterAll: H3K4me3_BRB --"))
corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3$H3K4me3, corr_H3p3clus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("H3K4me3", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
[1] "----- UI --- H3p3clusterAll: H3K4me3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 0h --- H3p3clusterAll: H3K4me3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 24h --- H3p3clusterAll: H3K4me3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 48h --- H3p3clusterAll: H3K4me3_BRB --"
タイのため正確な p 値を計算することができません
print("~~ H3K27ac_BRB ~~")
[1] "~~ H3K27ac_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clusterAll: H3K27ac_BRB --"))
corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3$H3K27ac, corr_H3p3clus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("H3K27ac", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
[1] "----- UI --- H3p3clusterAll: H3K27ac_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 0h --- H3p3clusterAll: H3K27ac_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 24h --- H3p3clusterAll: H3K27ac_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 48h --- H3p3clusterAll: H3K27ac_BRB --"
タイのため正確な p 値を計算することができません
print("~~ H3K27me3_BRB ~~")
[1] "~~ H3K27me3_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clusterAll: H3K27me3_BRB --"))
corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3$H3K27me3, corr_H3p3clus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("H3K27me3", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
[1] "----- UI --- H3p3clusterAll: H3K27me3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 0h --- H3p3clusterAll: H3K27me3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 24h --- H3p3clusterAll: H3K27me3_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 48h --- H3p3clusterAll: H3K27me3_BRB --"
タイのため正確な p 値を計算することができません
print("~~ ATAC_BRB ~~")
[1] "~~ ATAC_BRB ~~"
for (i in 1:length(Time_list)) {
if ((i == 1)|(i == 4)) {
print(paste("-----",Time_list[i], "--- H3p3clusterAll: ATAC_BRB --"))
corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3$ATAC, corr_H3p3clus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("ATAC", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
}
[1] "----- UI --- H3p3clusterAll: ATAC_BRB --"
タイのため正確な p 値を計算することができません
[1] "----- 48h --- H3p3clusterAll: ATAC_BRB --"
タイのため正確な p 値を計算することができません
print("~~ H3p3_ATAC ~~")
[1] "~~ H3p3_ATAC ~~"
for (i in 1:length(Time_list)) {
if ((i == 1)|(i == 4)) {
print(paste("-----",Time_list[i], "--- H3p3clusterAll: H3p3_ATAC --"))
corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3$H3p3, corr_H3p3clus3$ATAC, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3All",Compare=paste("H3p3", "ATAC",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
}
[1] "----- UI --- H3p3clusterAll: H3p3_ATAC --"
タイのため正確な p 値を計算することができません
[1] "----- 48h --- H3p3clusterAll: H3p3_ATAC --"
タイのため正確な p 値を計算することができません
cortest_result_s_H3p3clus3All <- cortest_result_s %>% group_by(target,time,Compare) %>% spread(key="Cor_test",value=Value) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% arrange(time)
cortest_result_s_H3p3clus3All <- cortest_result_s_H3p3clus3All %>% left_join(Count_FC_cutoff_H3p3clus3)
Joining, by = "time"
print(cortest_result_s_H3p3clus3All)
cortest_result_s_H3p3clus3All %>% readr::write_csv("./log2FC/tables/Cortest_result_spearman_H3p3clus3All.csv")
Time_list <- c("UI","0h","24h","48h")
Count_FC_cutoff_H3p3clus3_BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% f_gene_BRBclus3 %>% ungroup() %>% group_by(time) %>% summarise(Plot_genes=n()) #図中の数(BRB DEG cluster3のみ)
`summarise()` ungrouping output (override with `.groups` argument)
#######
print("~~ H3p3_BRB ~~")
[1] "~~ H3p3_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clus3BRBclus3: H3p3_BRB --"))
corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3BRBclus3$H3p3, corr_H3p3clus3BRBclus3$BRB, method="spearman")
#print(tttttt)
if (i == 1) {
cortest_result_s <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("H3p3", "BRB",sep="_"))
}
else {
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("H3p3", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
}
[1] "----- UI --- H3p3clus3BRBclus3: H3p3_BRB --"
[1] "----- 0h --- H3p3clus3BRBclus3: H3p3_BRB --"
[1] "----- 24h --- H3p3clus3BRBclus3: H3p3_BRB --"
[1] "----- 48h --- H3p3clus3BRBclus3: H3p3_BRB --"
print("~~ H3K4me3_BRB ~~")
[1] "~~ H3K4me3_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clus3BRBclus3: H3K4me3_BRB --"))
corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3BRBclus3$H3K4me3, corr_H3p3clus3BRBclus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("H3K4me3", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
[1] "----- UI --- H3p3clus3BRBclus3: H3K4me3_BRB --"
[1] "----- 0h --- H3p3clus3BRBclus3: H3K4me3_BRB --"
[1] "----- 24h --- H3p3clus3BRBclus3: H3K4me3_BRB --"
[1] "----- 48h --- H3p3clus3BRBclus3: H3K4me3_BRB --"
print("~~ H3K27ac_BRB ~~")
[1] "~~ H3K27ac_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clus3BRBclus3: H3K27ac_BRB --"))
corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3BRBclus3$H3K27ac, corr_H3p3clus3BRBclus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("H3K27ac", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
[1] "----- UI --- H3p3clus3BRBclus3: H3K27ac_BRB --"
[1] "----- 0h --- H3p3clus3BRBclus3: H3K27ac_BRB --"
[1] "----- 24h --- H3p3clus3BRBclus3: H3K27ac_BRB --"
[1] "----- 48h --- H3p3clus3BRBclus3: H3K27ac_BRB --"
print("~~ H3K27me3_BRB ~~")
[1] "~~ H3K27me3_BRB ~~"
for (i in 1:length(Time_list)) {
print(paste("-----",Time_list[i], "--- H3p3clus3BRBclus3: H3K27me3_BRB --"))
corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3BRBclus3$H3K27me3, corr_H3p3clus3BRBclus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("H3K27me3", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
[1] "----- UI --- H3p3clus3BRBclus3: H3K27me3_BRB --"
[1] "----- 0h --- H3p3clus3BRBclus3: H3K27me3_BRB --"
[1] "----- 24h --- H3p3clus3BRBclus3: H3K27me3_BRB --"
[1] "----- 48h --- H3p3clus3BRBclus3: H3K27me3_BRB --"
print("~~ ATAC_BRB ~~")
[1] "~~ ATAC_BRB ~~"
for (i in 1:length(Time_list)) {
if ((i == 1)|(i == 4)) {
print(paste("-----",Time_list[i], "--- H3p3clus3BRBclus3: ATAC_BRB --"))
corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3BRBclus3$ATAC, corr_H3p3clus3BRBclus3$BRB, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("ATAC", "BRB",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
}
[1] "----- UI --- H3p3clus3BRBclus3: ATAC_BRB --"
[1] "----- 48h --- H3p3clus3BRBclus3: ATAC_BRB --"
print("~~ H3p3_ATAC ~~")
[1] "~~ H3p3_ATAC ~~"
for (i in 1:length(Time_list)) {
if ((i == 1)|(i == 4)) {
print(paste("-----",Time_list[i], "--- H3p3clus3BRBclus3: H3p3_ATAC --"))
corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time==Time_list[i]))
tttttt <- cor.test(corr_H3p3clus3BRBclus3$H3p3, corr_H3p3clus3BRBclus3$ATAC, method="spearman")
#print(tttttt)
ssssss <- unlist(tttttt) %>% as.data.frame() %>% tibble::rownames_to_column("Cor_test") %>% as_tibble() %>% dplyr::rename(Value=".") %>% mutate(time=Time_list[i],target="H3p3clus3BRBclus3",Compare=paste("H3p3", "ATAC",sep="_"))
cortest_result_s <- bind_rows(cortest_result_s, ssssss)
}
}
[1] "----- UI --- H3p3clus3BRBclus3: H3p3_ATAC --"
[1] "----- 48h --- H3p3clus3BRBclus3: H3p3_ATAC --"
cortest_result_s_H3p3clus3BRBclus3 <- cortest_result_s %>% group_by(target,time,Compare) %>% spread(key="Cor_test",value=Value) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% arrange(time)
cortest_result_s_H3p3clus3BRBclus3 <- cortest_result_s_H3p3clus3BRBclus3 %>% left_join(Count_FC_cutoff_H3p3clus3_BRBclus3)
Joining, by = "time"
print(cortest_result_s_H3p3clus3BRBclus3)
cortest_result_s_H3p3clus3BRBclus3 %>% readr::write_csv("./log2FC/tables/Cortest_results_spearman_H3p3clus3BRBclus3.csv")
# All
rank_corr_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% group_by(time) %>% mutate(rank_H3p3=rank(H3p3), rank_H3K4me3=rank(H3K4me3), rank_H3K27ac=rank(H3K27ac), rank_H3K27me3=rank(H3K27me3), rank_ATAC=rank(ATAC), rank_BRB=rank(BRB))
#H3p3clus3BRBclus3
rank_corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3") %>% group_by(time) %>% mutate(rank_H3p3=rank(H3p3), rank_H3K4me3=rank(H3K4me3), rank_H3K27ac=rank(H3K27ac), rank_H3K27me3=rank(H3K27me3), rank_ATAC=rank(ATAC), rank_BRB=rank(BRB))
###
fcplot <- rank_corr_H3p3clus3BRBclus3 %>% ggplot(aes(y=rank_BRB, x=rank_H3p3)) + facet_wrap(~time,nrow=1, scales = "free") +
xlim(0, NA) + ylim(0, NA) + geom_point(alpha = 0.6, size=1.0)+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank())
#+ xlim(0,60) + ylim(0,60)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3p3_BRB"),aes(x=0,y=60,label=paste(sprintf("%4.3e", estimate.rho)," (",Plot_genes," genes)", sep="")), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 2.0)
#+ stat_smooth()
#+ stat_smooth(method = "lm", colour = "black", size = 1)
#+ geom_smooth(method = lm, se = FALSE)
fcplot
#rank_corr_H3p3clus3BRBclus3 <- plot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% filter((time=="48h")) %>% mutate(rank_H3p3=rank(H3p3), rank_H3K4me3=rank(H3K4me3), rank_H3K27ac=rank(H3K27ac), rank_H3K27me3=rank(H3K27me3), rank_ATAC=rank(ATAC), rank_BRB=rank(BRB))
#cor.test(rank_corr_H3p3clus3BRBclus3$rank_H3p3, rank_corr_H3p3clus3BRBclus3$rank_BRB)
### 並び替えのため
groups_BRB_arr <- c("BRB_UI_DoxMinus","BRB_UI_DoxPlus","BRB_0h_DoxMinus","BRB_0h_DoxPlus","BRB_24h_DoxMinus","BRB_24h_DoxPlus","BRB_48h_DoxMinus","BRB_48h_DoxPlus")
groups_ATAC_arr <- c("ATAC_UI_DoxMinus","ATAC_UI_DoxPlus","ATAC_48h_DoxMinus","ATAC_48h_DoxPlus")
groups_H3p3_arr <- c("H3p3_UI_DoxMinus","H3p3_UI_DoxPlus","H3p3_0h_DoxMinus","H3p3_0h_DoxPlus","H3p3_24h_DoxMinus","H3p3_24h_DoxPlus","H3p3_48h_DoxMinus","H3p3_48h_DoxPlus")
groups_H3K4me3_arr <- c("H3K4me3_UI_DoxMinus","H3K4me3_UI_DoxPlus","H3K4me3_0h_DoxMinus","H3K4me3_0h_DoxPlus","H3K4me3_24h_DoxMinus","H3K4me3_24h_DoxPlus","H3K4me3_48h_DoxMinus","H3K4me3_48h_DoxPlus")
groups_H3K27ac_arr <- c("H3K27ac_UI_DoxMinus","H3K27ac_UI_DoxPlus","H3K27ac_0h_DoxMinus","H3K27ac_0h_DoxPlus","H3K27ac_24h_DoxMinus","H3K27ac_24h_DoxPlus","H3K27ac_48h_DoxMinus","H3K27ac_48h_DoxPlus")
groups_H3K27me3_arr <- c("H3K27me3_UI_DoxMinus","H3K27me3_UI_DoxPlus","H3K27me3_0h_DoxMinus","H3K27me3_0h_DoxPlus","H3K27me3_24h_DoxMinus","H3K27me3_24h_DoxPlus","H3K27me3_48h_DoxMinus","H3K27me3_48h_DoxPlus")
groupt_BRB_arr <- c("BRB_DoxMinus","BRB_DoxPlus")
groupt_ATAC_arr <- c("ATAC_DoxMinus","ATAC_DoxPlus")
groupt_H3p3_arr <- c("H3p3_DoxMinus","H3p3_DoxPlus")
groupt_H3K4me3_arr <- c("H3K4me3_DoxMinus","H3K4me3_DoxPlus")
groupt_H3K27ac_arr <- c("H3K27ac_DoxMinus","H3K27ac_DoxPlus")
groupt_H3K27me3_arr <- c("H3K27me3_DoxMinus","H3K27me3_DoxPlus")
gggg_list2 <- c("ens_gene", "ext_gene", "time","BRBgroupMean", "H3p3cluster", "BRBDEGcluster", all_of(groupt_H3p3_arr),all_of(groupt_H3K4me3_arr),all_of(groupt_H3K27ac_arr),all_of(groupt_H3K27me3_arr),all_of(groupt_ATAC_arr),all_of(groupt_BRB_arr))
gggg_list1 <- c("ens_gene", "H3p3cluster", "BRBDEGcluster",all_of(groups_H3p3_arr),all_of(groups_H3K4me3_arr),all_of(groups_H3K27ac_arr),all_of(groups_H3K27me3_arr),all_of(groups_ATAC_arr),all_of(groups_BRB_arr))
## 全て結合
norm_chilatac_cutoff <- norm_plotlist_all %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene) %>% dplyr::select(ens_gene,sample,group,time,type,seq,rep,normalized)
norm_brb_cutoff <- norm_BRB %>% filter(ens_gene %in% norm_BRB_cutoff$ens_gene) %>% dplyr::select(ens_gene,sample,group,time,type,seq,rep,norm) %>% rename(normalized=norm)
## 全て結合
norm_plotlist_cutoff <- bind_rows(norm_chilatac_cutoff,norm_brb_cutoff) %>% mutate(seq=factor(seq, c("H3p3", "H3K4me3","H3K27ac","H3K27me3","ATAC","BRB"))) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) # ChIL & ATAC & BRB (cut off 全て)
norm_typemean_cutoff <- norm_plotlist_cutoff %>% group_by(ens_gene,group,time,type,seq) %>% summarise(TypeMean=mean(normalized)) %>% mutate(seq_type=paste(seq,type,sep="_"))
`summarise()` regrouping output by 'ens_gene', 'group', 'time', 'type' (override with `.groups` argument)
spread_norm_typemean_cutoff <- norm_typemean_cutoff %>% ungroup() %>% dplyr::select(ens_gene,time,seq_type,TypeMean) %>% mutate(seq_type=factor(seq_type, c("H3p3_DoxMinus","H3p3_DoxPlus", "H3K4me3_DoxMinus","H3K4me3_DoxPlus","H3K27ac_DoxMinus","H3K27ac_DoxPlus","H3K27me3_DoxMinus","H3K27me3_DoxPlus","ATAC_DoxMinus","ATAC_DoxPlus","BRB_DoxMinus","BRB_DoxPlus"))) %>% arrange(seq_type) %>% spread(key=seq_type,value=TypeMean)
spread_norm_typemean_cutoff <- spread_norm_typemean_cutoff %>% left_join(dplyr::select(norm_BRB_cutoff, ens_gene, ext_gene, time, groupMean) %>% dplyr::rename(BRBgroupMean=groupMean))
Joining, by = c("ens_gene", "time")
nrow(spread_norm_typemean_cutoff)
[1] 32636
spread_norm_typemean_cutoff_clus <- spread_norm_typemean_cutoff %>% left_join(dplyr::select(rrres_allH3p3,ens_gene,cluster)) %>% dplyr::rename(H3p3cluster=cluster) %>% left_join(dplyr::select(cluster_BRBlist,ens_gene,cluster)) %>% dplyr::rename(BRBDEGcluster=cluster) %>% dplyr::select(all_of(gggg_list2))
Joining, by = "ens_gene"
Joining, by = "ens_gene"
nrow(spread_norm_typemean_cutoff_clus)
[1] 32636
#####
time_norm_typemean_cutoff <- norm_typemean_cutoff %>% ungroup() %>% dplyr::select(ens_gene,group,TypeMean) %>% group_by(ens_gene) %>% spread(key=group,value=TypeMean)
nrow(time_norm_typemean_cutoff)
[1] 8159
time_norm_typemean_cutoff_clus <- time_norm_typemean_cutoff %>% left_join(dplyr::select(rrres_allH3p3,ens_gene,cluster)) %>% dplyr::rename(H3p3cluster=cluster) %>% left_join(dplyr::select(cluster_BRBlist,ens_gene,cluster)) %>% dplyr::rename(BRBDEGcluster=cluster) %>% dplyr::select(all_of(gggg_list1))
Joining, by = "ens_gene"
Joining, by = "ens_gene"
nrow(time_norm_typemean_cutoff_clus)
[1] 8159
#####
filename <- "./Correlation/tables/normTypeMean_All_cutoff10.csv"
print(filename)
[1] "./Correlation/tables/normTypeMean_All_cutoff10.csv"
readr::write_csv(spread_norm_typemean_cutoff_clus,filename)
head(spread_norm_typemean_cutoff_clus)
nrow(spread_norm_typemean_cutoff_clus)
[1] 32636
filename <- "./Correlation/tables/normTypeMean_All_cutoff10__timever.csv"
print(filename)
[1] "./Correlation/tables/normTypeMean_All_cutoff10__timever.csv"
readr::write_csv(time_norm_typemean_cutoff_clus,filename)
head(time_norm_typemean_cutoff_clus)
nrow(time_norm_typemean_cutoff_clus)
[1] 8159
#f_gene_H3p3clus3 <- function(x) x %>% filter(H3p3cluster=="3")
#f_gene_BRBclus3 <- function(x) x %>% filter(BRBDEGcluster=="3")
#list_gene_qpcr <- c("Acta1","Myh3","Ttn","Myog")
corr_typemean_cutoff_H3p3clus3 <- time_norm_typemean_cutoff_clus %>% ungroup() %>% f_gene_H3p3clus3
nrow(corr_typemean_cutoff_H3p3clus3)
[1] 3363
corr_typemean_cutoff_H3p3clus3BRBclus3 <- time_norm_typemean_cutoff_clus %>% ungroup() %>% f_gene_H3p3clus3 %>% f_gene_BRBclus3
nrow(corr_typemean_cutoff_H3p3clus3BRBclus3)
[1] 55
####
corr_typemean_cutoff_H3p3clus3_filter <- spread_norm_typemean_cutoff_clus %>% filter(!is.na(BRBgroupMean)) %>% ungroup() %>% f_gene_H3p3clus3
nrow(corr_typemean_cutoff_H3p3clus3_filter)
[1] 11802
corr_typemean_cutoff_H3p3clus3BRBclus3_filter <- spread_norm_typemean_cutoff_clus %>% filter(!is.na(BRBgroupMean)) %>% ungroup() %>% f_gene_H3p3clus3 %>% f_gene_BRBclus3
nrow(corr_typemean_cutoff_H3p3clus3BRBclus3_filter)
[1] 190
corr_typemean_cutoff_H3p3clus3_filter %>% group_by(time) %>% summarise(count=n())
`summarise()` ungrouping output (override with `.groups` argument)
corr_typemean_cutoff_H3p3clus3BRBclus3_filter %>% group_by(time) %>% summarise(count=n())
`summarise()` ungrouping output (override with `.groups` argument)
#filename <- "./log2FC/tables/Spread_log2FC_ChILATACBRB_cutoff10__withCluster__H3p3clus3.csv"
#print(filename)
#readr::write_csv(spread_norm_typemean_cutoff_H3p3clus3,filename)
#spread_norm_typemean_cutoff_H3p3clus3
#spread_norm_typemean_cutoff_H3p3clus3 %>% f_gene_H3p3clus3 %>% f_gene_BRBclus3
#spread_norm_typemean_cutoff_H3p3clus3 %>% f_gene_H3p3clus3 %>% f_gene_BRBclus3 %>% filter(ext_gene %in% list_gene_qpcr)
#H3p3clus3cutoff <- spread_norm_typemean_cutoff_H3p3clus3 %>% ungroup() %>% dplyr::select(ens_gene) %>% unique() %>% nrow()
#H3p3clus3cutoff_brbclus3 <- spread_norm_typemean_cutoff_H3p3clus3 %>% f_gene_BRBclus3 %>% ungroup() %>% dplyr::select(ens_gene) %>% unique() %>% nrow()
#nrow(z_H3p3clus3)
#H3p3clus3cutoff
#H3p3clus3cutoff_brbclus3
いずれかで BRB norm (group) > 10を満たすものでプロット
library(corrplot)
#### calculate correlation
print("Genes list")
[1] "Genes list"
nrow(corr_typemean_cutoff_H3p3clus3)
[1] 3363
nrow(corr_typemean_cutoff_H3p3clus3BRBclus3)
[1] 55
mydata.cor.All.alltime <- corr_typemean_cutoff_H3p3clus3 %>% ungroup() %>% dplyr::select(-"ens_gene", -"H3p3cluster", -"BRBDEGcluster") %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.alltime <- corr_typemean_cutoff_H3p3clus3BRBclus3 %>% ungroup() %>% dplyr::select(-"ens_gene", -"H3p3cluster", -"BRBDEGcluster") %>% cor(method = c("spearman"))
cor_All_alltime <- as.data.frame(mydata.cor.All.alltime) %>% tibble::rownames_to_column("group") %>% as_tibble
cor_BRBclus3_alltime <- as.data.frame(mydata.cor.BRBclus3.alltime) %>% tibble::rownames_to_column("group") %>% as_tibble
####
title_1 <- paste("H3.3 Cluster3 & BRB DEG Cluster3:",nrow(corr_typemean_cutoff_H3p3clus3BRBclus3),"genes",sep=" ")
title_2 <- paste("H3.3 Cluster3:",nrow(corr_typemean_cutoff_H3p3clus3),"genes",sep=" ")
breaksList = seq(-1, 1, by = 0.05)
Color__a0 <- rev(brewer.pal(n = 11, name = "RdYlBu"))
Color__a <- colorRampPalette(Color__a0)(length(breaksList))
gaps_1 <- c(8,16,24,32,36)
pheatmap::pheatmap(mydata.cor.BRBclus3.alltime, main = title_1, scale = "none", color = Color__a,breaks = breaksList, border_color=NA)
pheatmap::pheatmap(mydata.cor.All.alltime, main = title_2, scale = "none", color = Color__a,breaks = breaksList, border_color=NA)
pheatmap::pheatmap(mydata.cor.BRBclus3.alltime, main = title_1, scale = "none", color = Color__a,breaks = breaksList, border_color=NA, cluster_rows = FALSE, cluster_cols = FALSE, gaps_col=gaps_1, gaps_row=gaps_1)
pheatmap::pheatmap(mydata.cor.All.alltime, main = title_2, scale = "none", color = Color__a,breaks = breaksList, border_color=NA, cluster_rows = FALSE, cluster_cols = FALSE, gaps_col=gaps_1, gaps_row=gaps_1)
corrplot(mydata.cor.BRBclus3.alltime, diag = FALSE, col = Color__a)
corrplot(mydata.cor.All.alltime, diag = FALSE, col = Color__a)
#####
filename <- "./Correlation/tables/Cortest_normTypeMean_spearman_H3p3clus3All_cutoff10__timever.csv"
print(filename)
[1] "./Correlation/tables/Cortest_normTypeMean_spearman_H3p3clus3All_cutoff10__timever.csv"
readr::write_csv(cor_All_alltime,filename)
print(cor_All_alltime)
nrow(cor_All_alltime)
[1] 44
filename <- "./Correlation/tables/Cortest_normTypeMean_spearman_H3p3clus3BRBclus3_cutoff10__timever.csv"
print(filename)
[1] "./Correlation/tables/Cortest_normTypeMean_spearman_H3p3clus3BRBclus3_cutoff10__timever.csv"
readr::write_csv(cor_BRBclus3_alltime,filename)
print(cor_BRBclus3_alltime)
nrow(cor_BRBclus3_alltime)
[1] 44
#corrplot(mydata.cor.BRBclus3.alltime, diag = FALSE, type = "upper", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.alltime, diag = FALSE, type = "lower", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.alltime, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.alltime, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
## 表示リスト
gggg_plot_corlist <- c(all_of(groups_H3p3_arr),all_of(groups_H3K4me3_arr),all_of(groups_H3K27ac_arr),all_of(groups_BRB_arr))
cor_All_alltime_select <- cor_All_alltime %>% dplyr::select(group, all_of(gggg_plot_corlist)) %>% filter(group %in% gggg_plot_corlist)
mat_cor_All_alltime_select <- cor_All_alltime_select %>% dplyr::select(-group) %>% as.matrix()
rownames(mat_cor_All_alltime_select) <- cor_All_alltime_select$group
cor_BRBclus3_alltime_select <- cor_BRBclus3_alltime %>% dplyr::select(group, all_of(gggg_plot_corlist)) %>% filter(group %in% gggg_plot_corlist)
mat_cor_BRBclus3_alltime_select <- cor_BRBclus3_alltime_select %>% dplyr::select(-group) %>% as.matrix()
rownames(mat_cor_BRBclus3_alltime_select) <- cor_BRBclus3_alltime_select$group
gaps_2 <- c(8,16,24)
#### plot (select only)
#breaksList = seq(-1, 1, by = 0.05)
#Color__a0 <- rev(brewer.pal(n = 11, name = "RdYlBu"))
#Color__a <- colorRampPalette(Color__a0)(length(breaksList))
pheatmap::pheatmap(mat_cor_BRBclus3_alltime_select, main = title_1, scale = "none", color = Color__a,breaks = breaksList, border_color=NA)
pheatmap::pheatmap(mat_cor_All_alltime_select, main = title_2, scale = "none", color = Color__a,breaks = breaksList, border_color=NA)
pheatmap::pheatmap(mat_cor_BRBclus3_alltime_select, main = title_1, scale = "none", color = Color__a,breaks = breaksList, border_color=NA, cluster_rows = FALSE, cluster_cols = FALSE, gaps_col=gaps_2, gaps_row=gaps_2)
pheatmap::pheatmap(mat_cor_All_alltime_select, main = title_2, scale = "none", color = Color__a,breaks = breaksList, border_color=NA, cluster_rows = FALSE, cluster_cols = FALSE, gaps_col=gaps_2, gaps_row=gaps_2)
corrplot(mat_cor_BRBclus3_alltime_select, diag = FALSE, col = Color__a)
corrplot(mat_cor_All_alltime_select, diag = FALSE, col = Color__a)
groups_UI_arr <- c(
"H3p3_UI_DoxMinus","H3p3_UI_DoxPlus",
"H3K4me3_UI_DoxMinus","H3K4me3_UI_DoxPlus",
"H3K27ac_UI_DoxMinus","H3K27ac_UI_DoxPlus",
"H3K27me3_UI_DoxMinus","H3K27me3_UI_DoxPlus",
"ATAC_UI_DoxMinus","ATAC_UI_DoxPlus",
"BRB_UI_DoxMinus","BRB_UI_DoxPlus")
groups_0h_arr <- c(
"H3p3_0h_DoxMinus","H3p3_0h_DoxPlus",
"H3K4me3_0h_DoxMinus","H3K4me3_0h_DoxPlus",
"H3K27ac_0h_DoxMinus","H3K27ac_0h_DoxPlus",
"H3K27me3_0h_DoxMinus","H3K27me3_0h_DoxPlus",
"BRB_0h_DoxMinus","BRB_0h_DoxPlus")
groups_24h_arr <- c(
"H3p3_24h_DoxMinus","H3p3_24h_DoxPlus",
"H3K4me3_24h_DoxMinus","H3K4me3_24h_DoxPlus",
"H3K27ac_24h_DoxMinus","H3K27ac_24h_DoxPlus",
"H3K27me3_24h_DoxMinus","H3K27me3_24h_DoxPlus",
"BRB_24h_DoxMinus","BRB_24h_DoxPlus")
groups_48h_arr <- c(
"H3p3_48h_DoxMinus","H3p3_48h_DoxPlus",
"H3K4me3_48h_DoxMinus","H3K4me3_48h_DoxPlus",
"H3K27ac_48h_DoxMinus","H3K27ac_48h_DoxPlus",
"H3K27me3_48h_DoxMinus","H3K27me3_48h_DoxPlus",
"ATAC_48h_DoxMinus","ATAC_48h_DoxPlus",
"BRB_48h_DoxMinus","BRB_48h_DoxPlus")
groups_UI_arr_s <- c(
"H3p3_UI_DoxMinus","H3p3_UI_DoxPlus",
"H3K4me3_UI_DoxMinus","H3K4me3_UI_DoxPlus",
"H3K27ac_UI_DoxMinus","H3K27ac_UI_DoxPlus",
"BRB_UI_DoxMinus","BRB_UI_DoxPlus")
groups_0h_arr_s <- c(
"H3p3_0h_DoxMinus","H3p3_0h_DoxPlus",
"H3K4me3_0h_DoxMinus","H3K4me3_0h_DoxPlus",
"H3K27ac_0h_DoxMinus","H3K27ac_0h_DoxPlus",
"BRB_0h_DoxMinus","BRB_0h_DoxPlus")
groups_24h_arr_s <- c(
"H3p3_24h_DoxMinus","H3p3_24h_DoxPlus",
"H3K4me3_24h_DoxMinus","H3K4me3_24h_DoxPlus",
"H3K27ac_24h_DoxMinus","H3K27ac_24h_DoxPlus",
"BRB_24h_DoxMinus","BRB_24h_DoxPlus")
groups_48h_arr_s <- c(
"H3p3_48h_DoxMinus","H3p3_48h_DoxPlus",
"H3K4me3_48h_DoxMinus","H3K4me3_48h_DoxPlus",
"H3K27ac_48h_DoxMinus","H3K27ac_48h_DoxPlus",
"BRB_48h_DoxMinus","BRB_48h_DoxPlus")
## 表示リスト
gggg_plot_corlist_time <- c(all_of(groups_UI_arr_s),all_of(groups_0h_arr_s),all_of(groups_24h_arr_s),all_of(groups_48h_arr_s))
cor_All_alltime_select2 <- cor_All_alltime %>% dplyr::select(group, all_of(gggg_plot_corlist_time)) %>% filter(group %in% gggg_plot_corlist_time) %>% mutate(group=factor(group,gggg_plot_corlist_time)) %>% arrange(group)
mat_cor_All_alltime_select2 <- cor_All_alltime_select2 %>% dplyr::select(-group) %>% as.matrix()
rownames(mat_cor_All_alltime_select2) <- cor_All_alltime_select2$group
cor_BRBclus3_alltime_select2 <- cor_BRBclus3_alltime%>% dplyr::select(group, all_of(gggg_plot_corlist_time)) %>% filter(group %in% gggg_plot_corlist_time) %>% mutate(group=factor(group,gggg_plot_corlist_time)) %>% arrange(group)
mat_cor_BRBclus3_alltime_select2 <- cor_BRBclus3_alltime_select2 %>% dplyr::select(-group) %>% as.matrix()
rownames(mat_cor_BRBclus3_alltime_select2) <- cor_BRBclus3_alltime_select2$group
gaps_2 <- c(8,16,24)
#### plot (select only)
pheatmap::pheatmap(mat_cor_BRBclus3_alltime_select2, main = title_1, scale = "none", color = Color__a,breaks = breaksList, border_color=NA, cluster_rows = FALSE, cluster_cols = FALSE, gaps_col=gaps_2, gaps_row=gaps_2)
pheatmap::pheatmap(mat_cor_All_alltime_select2, main = title_2, scale = "none", color = Color__a,breaks = breaksList, border_color=NA, cluster_rows = FALSE, cluster_cols = FALSE, gaps_col=gaps_2, gaps_row=gaps_2)
NA
NA
NA
mydata.cor.BRBclus3.0h <- corr_typemean_cutoff_H3p3clus3BRBclus3_filter %>% ungroup() %>% filter(time=="0h") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.24h <- corr_typemean_cutoff_H3p3clus3BRBclus3_filter %>% ungroup() %>% filter(time=="24h") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.48h <- corr_typemean_cutoff_H3p3clus3BRBclus3_filter %>% ungroup() %>% filter(time=="48h") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
mydata.cor.All.UI <- corr_typemean_cutoff_H3p3clus3_filter %>% ungroup() %>% filter(time=="UI") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
mydata.cor.All.0h <- corr_typemean_cutoff_H3p3clus3_filter %>% ungroup() %>% filter(time=="0h") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
mydata.cor.All.24h <- corr_typemean_cutoff_H3p3clus3_filter %>% ungroup() %>% filter(time=="24h") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
mydata.cor.All.48h <- corr_typemean_cutoff_H3p3clus3_filter %>% ungroup() %>% filter(time=="48h") %>% dplyr::select(BRB_DoxMinus, BRB_DoxPlus, H3p3_DoxMinus,H3p3_DoxPlus,H3K4me3_DoxMinus,H3K4me3_DoxPlus,H3K27ac_DoxMinus,H3K27ac_DoxPlus) %>% cor(method = c("spearman"))
#corrplot(mydata.cor.BRBclus3.UI, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.0h, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.24h, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.48h, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.UI, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.0h, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.24h, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.48h, diag = FALSE, method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.BRBclus3.UI, diag = FALSE, method = "color")
corrplot(mydata.cor.BRBclus3.0h, diag = FALSE, method = "color")
corrplot(mydata.cor.BRBclus3.24h, diag = FALSE, method = "color")
corrplot(mydata.cor.BRBclus3.48h, diag = FALSE, method = "color")
corrplot(mydata.cor.All.UI, diag = FALSE, method = "color")
corrplot(mydata.cor.All.0h, diag = FALSE, method = "color")
corrplot(mydata.cor.All.24h, diag = FALSE, method = "color")
corrplot(mydata.cor.All.48h, diag = FALSE, method = "color")
corrplot(mydata.cor.BRBclus3.UI, diag = FALSE)
corrplot(mydata.cor.BRBclus3.0h, diag = FALSE)
corrplot(mydata.cor.BRBclus3.24h, diag = FALSE)
corrplot(mydata.cor.BRBclus3.48h, diag = FALSE)
corrplot(mydata.cor.All.UI, diag = FALSE)
corrplot(mydata.cor.All.0h, diag = FALSE)
corrplot(mydata.cor.All.24h, diag = FALSE)
corrplot(mydata.cor.All.48h, diag = FALSE)
#corrplot(mydata.cor.BRBclus3.UI, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.BRBclus3.0h, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.BRBclus3.24h, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.BRBclus3.48h, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.All.UI, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.All.0h, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.All.24h, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot(mydata.cor.All.48h, diag = FALSE, method = "color", col = cm.colors(100))
#corrplot.mixed(cor(mydata.cor), order="hclust", tl.col="black")
#pheatmap::pheatmap(mydata.cor.All.UI,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.All.0h,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.All.24h,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.All.48h,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.BRBclus3.UI,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.BRBclus3.0h,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.BRBclus3.24h,color=viridis::viridis(256),scale = "none")
#pheatmap::pheatmap(mydata.cor.BRBclus3.48h,color=viridis::viridis(256),scale = "none")
#corrplot(mydata.cor.BRBclus3.UI, diag = FALSE, type = "upper", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.0h, diag = FALSE, type = "upper", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.24h, diag = FALSE, type = "upper", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.48h, diag = FALSE, type = "upper", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.UI, diag = FALSE, type = "lower", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.0h, diag = FALSE, type = "lower", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.24h, diag = FALSE, type = "lower", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.48h, diag = FALSE, type = "lower", method = "color", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.UI, diag = FALSE, type = "upper", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.0h, diag = FALSE, type = "upper", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.24h, diag = FALSE, type = "upper", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.BRBclus3.48h, diag = FALSE, type = "upper", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.UI, diag = FALSE, type = "lower", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.0h, diag = FALSE, type = "lower", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.24h, diag = FALSE, type = "lower", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot(mydata.cor.All.48h, diag = FALSE, type = "lower", method = "square", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot.mixed(mydata.cor.BRBclus3.UI, lower.col = "black")
#corrplot.mixed(mydata.cor.BRBclus3.0h, lower.col = "black")
#corrplot.mixed(mydata.cor.BRBclus3.24h, lower.col = "black")
#corrplot.mixed(mydata.cor.BRBclus3.48h, lower.col = "black")
#corrplot(mydata.cor,palette = "PuOr")
library(corrplot)
# = read.csv("https://wiki.q-researchsoftware.com/images/b/b9/Ownership.csv", header = TRUE, fileEncoding="latin1")
#mydata.cor = cor(mydata)
#mydata.cor = cor(mydata, method = c("spearman")
#corrplot(mydata.cor)
#mydata.cor.UI <- rank_corr_H3p3clus3BRBclus3 %>% filter(time=="UI") %>% dplyr::select(ens_gene,ext_gene,time,H3p3,H3K4me3,H3K27ac,BRB) %>% tidyr::gather(key="seq",value=FC,-ens_gene,-ext_gene,-time) %>% mutate(time_seq=paste(time,seq,sep="_")) %>% ungroup %>% dplyr::select(ens_gene,time_seq,FC) %>% tidyr::spread(key = time_seq,value=FC) %>% dplyr::select(-ens_gene) %>% cor(method = c("spearman"))
#mydata.cor.0h <- rank_corr_H3p3clus3BRBclus3 %>% filter(time=="0h") %>% dplyr::select(ens_gene,ext_gene,time,H3p3,H3K4me3,H3K27ac,BRB) %>% tidyr::gather(key="seq",value=FC,-ens_gene,-ext_gene,-time) %>% mutate(time_seq=paste(time,seq,sep="_")) %>% ungroup %>% dplyr::select(ens_gene,time_seq,FC) %>% tidyr::spread(key = time_seq,value=FC) %>% dplyr::select(-ens_gene) %>% cor(method = c("spearman"))
#mydata.cor.24h <- rank_corr_H3p3clus3BRBclus3 %>% filter(time=="24h") %>% dplyr::select(ens_gene,ext_gene,time,H3p3,H3K4me3,H3K27ac,BRB) %>% tidyr::gather(key="seq",value=FC,-ens_gene,-ext_gene,-time) %>% mutate(time_seq=paste(time,seq,sep="_")) %>% ungroup %>% dplyr::select(ens_gene,time_seq,FC) %>% tidyr::spread(key = time_seq,value=FC) %>% dplyr::select(-ens_gene) %>% cor(method = c("spearman"))
#mydata.cor.48h <- rank_corr_H3p3clus3BRBclus3 %>% filter(time=="48h") %>% dplyr::select(ens_gene,ext_gene,time,H3p3,H3K4me3,H3K27ac,BRB) %>% tidyr::gather(key="seq",value=FC,-ens_gene,-ext_gene,-time) %>% mutate(time_seq=paste(time,seq,sep="_")) %>% ungroup %>% dplyr::select(ens_gene,time_seq,FC) %>% tidyr::spread(key = time_seq,value=FC) %>% dplyr::select(-ens_gene) %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.UI <- rank_corr_H3p3clus3BRBclus3 %>% ungroup() %>% filter(time=="UI") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.0h <- rank_corr_H3p3clus3BRBclus3 %>% ungroup() %>% filter(time=="0h") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.24h <- rank_corr_H3p3clus3BRBclus3 %>% ungroup() %>% filter(time=="24h") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.BRBclus3.48h <- rank_corr_H3p3clus3BRBclus3 %>% ungroup() %>% filter(time=="48h") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.All.UI <- rank_corr_H3p3clus3 %>% ungroup() %>% filter(time=="UI") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.All.0h <- rank_corr_H3p3clus3 %>% ungroup() %>% filter(time=="0h") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.All.24h <- rank_corr_H3p3clus3 %>% ungroup() %>% filter(time=="24h") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
mydata.cor.All.48h <- rank_corr_H3p3clus3 %>% ungroup() %>% filter(time=="48h") %>% dplyr::select(BRB, H3p3,H3K4me3,H3K27ac) %>% cor(method = c("spearman"))
#corrplot.mixed(cor(mydata.cor), order="hclust", tl.col="black")
corrplot(mydata.cor.BRBclus3.UI, diag = FALSE, type = "upper", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.BRBclus3.0h, diag = FALSE, type = "upper", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.BRBclus3.24h, diag = FALSE, type = "upper", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.BRBclus3.48h, diag = FALSE, type = "upper", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.All.UI, diag = FALSE, type = "lower", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.All.0h, diag = FALSE, type = "lower", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.All.24h, diag = FALSE, type = "lower", col = rev(brewer.pal(n = 10, name = "RdBu")))
corrplot(mydata.cor.All.48h, diag = FALSE, type = "lower", col = rev(brewer.pal(n = 10, name = "RdBu")))
#corrplot.mixed(mydata.cor.BRBclus3.UI, lower.col = "black")
#corrplot.mixed(mydata.cor.BRBclus3.0h, lower.col = "black")
#corrplot.mixed(mydata.cor.BRBclus3.24h, lower.col = "black")
#corrplot.mixed(mydata.cor.BRBclus3.48h, lower.col = "black")
#corrplot(mydata.cor,palette = "PuOr")
spearmanとranking=>peason の相関係数が合っているかを調べること
#%>% group_by(target,time,Compare)
# group_by(aspect,gs_cat,gs_subcat) %>%
# mutate(padj=p.adjust(pval,"BH")) %>% ungroup()
#Cortest_H3p3clus3All <- readr::read_csv("./log2FC/tables/Cortest_result_spearman_H3p3clus3All.csv") %>% mutate(text=paste(" All (",Plot_genes,") Cor: ",sprintf("%4.3e", estimate.cor),", p.val: ",sprintf("%4.3e", p.value),sep="")) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% arrange(time)
#Cortest_H3p3clus3BRBclus3 <- readr::read_csv("./log2FC/tables/Cortest_result_spearman_H3p3clus3BRBclus3.csv") %>% mutate(text=paste("Clus3 (",Plot_genes,") Cor: ",sprintf("%4.3e", estimate.cor),", p.val: ",sprintf("%4.3e", p.value),sep="")) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% arrange(time)
Cortest_H3p3clus3All <- readr::read_csv("./log2FC/tables/Cortest_result_spearman_H3p3clus3All.csv") %>% mutate(text=paste(" All (",Plot_genes,") Spearman Cor: ",sprintf("%4.3e", estimate.rho),", p.val: ",sprintf("%4.3e", p.value),sep="")) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% arrange(time)
Cortest_H3p3clus3BRBclus3 <- readr::read_csv("./log2FC/tables/Cortest_results_spearman_H3p3clus3BRBclus3.csv") %>% mutate(text=paste("Clus3 (",Plot_genes,") Spearman Cor: ",sprintf("%4.3e", estimate.rho),", p.val: ",sprintf("%4.3e", p.value),sep="")) %>% mutate(time=factor(time, c("UI", "0h","24h","48h"))) %>% arrange(time)
Cortest_H3p3clus3All %>% readr::write_csv("./log2FC/tables/Cortest_result_spearman_H3p3clus3All_forPlot.csv")
Cortest_H3p3clus3BRBclus3 %>% readr::write_csv("./log2FC/tables/Cortest_result_spearman_H3p3clus3BRBclus3_forPlot.csv")
Showplot_all_FC_cutoff_H3p3clus3 <- plot_all_FC_cutoff_H3p3clus3 %>% mutate(H3p3vsBRB=case_when(((abs(BRB)>2)|(abs(H3p3)>1.0)|is.na(H3p3))~"Not Shown",TRUE~"Shown"),H3K4me3vsBRB=case_when(((abs(BRB)>2)|(abs(H3K4me3)>1.5)|is.na(H3K4me3))~"Not Shown",TRUE~"Shown"),H3K27acvsBRB=case_when(((abs(BRB)>2)|(abs(H3K27ac)>1.0)|is.na(H3K27ac))~"Not Shown",TRUE~"Shown"),H3K27me3vsBRB=case_when(((abs(BRB)>2)|(abs(H3K27me3)>1.0)|is.na(H3K27me3))~"Not Shown",TRUE~"Shown"),ATACvsBRB=case_when(((abs(BRB)>2)|(abs(ATAC)>0.4)|is.na(ATAC))~"Not Shown",TRUE~"Shown"))
Show_H3p3vsBRB_H3p3clus3All <- Showplot_all_FC_cutoff_H3p3clus3 %>% ungroup %>% group_by(H3p3vsBRB,time) %>% summarize(count=n()) %>% rename(Plot=H3p3vsBRB) %>% mutate(Compare="H3p3_BRB")
Show_H3K4me3vsBRB_H3p3clus3All <- Showplot_all_FC_cutoff_H3p3clus3 %>% ungroup %>% group_by(H3K4me3vsBRB,time) %>% summarize(count=n()) %>% rename(Plot=H3K4me3vsBRB) %>% mutate(Compare="H3K4me3_BRB")
Show_H3K27acvsBRB_H3p3clus3All <- Showplot_all_FC_cutoff_H3p3clus3 %>% ungroup %>% group_by(H3K27acvsBRB,time) %>% summarize(count=n()) %>% rename(Plot=H3K27acvsBRB) %>% mutate(Compare="H3K27ac_BRB")
Show_H3K27me3vsBRB_H3p3clus3All <- Showplot_all_FC_cutoff_H3p3clus3 %>% ungroup %>% group_by(H3K27me3vsBRB,time) %>% summarize(count=n()) %>% rename(Plot=H3K27me3vsBRB)%>% mutate(Compare="H3K27me3_BRB")
Show_ATACvsBRB_H3p3clus3All <- Showplot_all_FC_cutoff_H3p3clus3 %>% ungroup %>% group_by(ATACvsBRB,time) %>% summarize(count=n()) %>% rename(Plot=ATACvsBRB)%>% mutate(Compare="ATAC_BRB")
Show__vsBRB_H3p3clus3All <- bind_rows(Show_H3p3vsBRB_H3p3clus3All,Show_H3K4me3vsBRB_H3p3clus3All) %>% bind_rows(Show_H3K27acvsBRB_H3p3clus3All) %>% bind_rows(Show_H3K27me3vsBRB_H3p3clus3All) %>% bind_rows(Show_ATACvsBRB_H3p3clus3All) %>% mutate(target="H3p3clus3All") %>% mutate(Count=case_when(Plot=="Not Shown"~paste("(",count,")",sep=""),TRUE~as.character(count))) %>% mutate(Plot=factor(Plot, c("Shown","Not Shown"))) %>% arrange(Compare,time,Plot)
Show_H3p3vsBRB_H3p3clus3BRBclus3 <- Showplot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3") %>% ungroup %>% group_by(H3p3vsBRB,time) %>% summarize(count=n(),gene=paste(ext_gene,collapse = ",")) %>% rename(Plot=H3p3vsBRB) %>% mutate(Compare="H3p3_BRB")
Show_H3K4me3vsBRB_H3p3clus3BRBclus3 <- Showplot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% ungroup %>% group_by(H3K4me3vsBRB,time) %>% summarize(count=n(),gene=paste(ext_gene,collapse = ",")) %>% rename(Plot=H3K4me3vsBRB) %>% mutate(Compare="H3K4me3_BRB")
Show_H3K27acvsBRB_H3p3clus3BRBclus3 <- Showplot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% ungroup %>% group_by(H3K27acvsBRB,time) %>% summarize(count=n(),gene=paste(ext_gene,collapse = ",")) %>% rename(Plot=H3K27acvsBRB) %>% mutate(Compare="H3K27ac_BRB")
Show_H3K27me3vsBRB_H3p3clus3BRBclus3 <- Showplot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% ungroup %>% group_by(H3K27me3vsBRB,time) %>% summarize(count=n(),gene=paste(ext_gene,collapse = ",")) %>% rename(Plot=H3K27me3vsBRB) %>% mutate(Compare="H3K27me3_BRB")
Show_ATACvsBRB_H3p3clus3BRBclus3 <- Showplot_all_FC_cutoff_H3p3clus3 %>% filter(BRBDEGcluster=="3")%>% ungroup %>% group_by(ATACvsBRB,time) %>% summarize(count=n(),gene=paste(ext_gene,collapse = ",")) %>% rename(Plot=ATACvsBRB) %>% mutate(Compare="ATAC_BRB")
Show__vsBRB_H3p3clus3BRBclus3 <- bind_rows(Show_H3p3vsBRB_H3p3clus3BRBclus3,Show_H3K4me3vsBRB_H3p3clus3BRBclus3) %>% bind_rows(Show_H3K27acvsBRB_H3p3clus3BRBclus3) %>% bind_rows(Show_H3K27me3vsBRB_H3p3clus3BRBclus3) %>% bind_rows(Show_ATACvsBRB_H3p3clus3BRBclus3) %>% mutate(target="H3p3clus3BRBclus3") %>% mutate(Count=case_when(Plot=="Not Shown"~paste("(",count,")",sep=""),TRUE~as.character(count))) %>% mutate(Plot=factor(Plot, c("Shown","Not Shown"))) %>% arrange(Compare,time,Plot)
####
Show__vsBRB_H3p3clus3All
Show__vsBRB_H3p3clus3BRBclus3
Show__vsBRB_H3p3clus3BRBclus3 %>% filter(Plot=="Not Shown")
Show__vsBRB_H3p3clus3All %>% readr::write_csv("./log2FC/tables/PLotshow_H3p3clus3All.csv")
Show__vsBRB_H3p3clus3BRBclus3 %>% readr::write_csv("./log2FC/tables/PLotshow_H3p3clus3BRBclus3.csv")
summa__vsBRB_H3p3clus3All <- Show__vsBRB_H3p3clus3All %>% ungroup() %>% group_by(target, time, Compare) %>% summarize(Show=paste(Count,collapse = " "))
summa__vsBRB_H3p3clus3BRBclus3 <- Show__vsBRB_H3p3clus3BRBclus3 %>% ungroup() %>% group_by(target, time, Compare) %>% summarize(Show=paste(Count,collapse = " "))
summa__vsBRB_H3p3clus3All %>% readr::write_csv("./log2FC/tables/PLotshow_H3p3clus3All_summary.csv")
summa__vsBRB_H3p3clus3BRBclus3 %>% readr::write_csv("./log2FC/tables/PLotshow_H3p3clus3BRBclus3_summary.csv")
summa__vsBRB_H3p3clus3_ALL_BRBclus3 <- bind_rows(summa__vsBRB_H3p3clus3All,summa__vsBRB_H3p3clus3BRBclus3) %>% mutate(target=factor(target, c("H3p3clus3All","H3p3clus3BRBclus3"))) %>% arrange(Compare,time,target) %>% mutate(clus=gsub("H3p3clus3","",target)) %>% ungroup() %>% group_by(time,Compare) %>% summarize(show=paste(Show,collapse = " / "),Cluster=paste(clus,collapse = " / "))
summa__vsBRB_H3p3clus3_ALL_BRBclus3 %>% readr::write_csv("./log2FC/tables/PLotshow_H3p3clus3_All_and_BRBclus3_summary.csv")
summa__vsBRB_H3p3clus3_ALL_BRBclus3
#density_color_low <- "#ECE038"
density_color_low <- "#FFFFFF"
#density_color_low <- "#ECE038"
density_color_low <- "#FFFFFF"
#density_color_high <- "#377EB8"
density_color_high <- "blue"
#density_color_low <- #FFFFFF"
#density_color_mid <- "yellow"
#density_color_high <- "red"
binsize <- 7
pppplottitle <- paste("log2 FC (Dox + vs -)\nBRB normalized count (Time, avg) > ",Set_cutoff,"\n H3.3 clus3: ",nrow(z_H3p3clus3)," genes\n Plot: ",H3p3clus3cutoff," genes\n BRB clus3: ",H3p3clus3cutoff_brbclus3," genes",sep="")
###
fcplot <-plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3p3)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + stat_ellipse(type = "norm",color=density_color_high,size=0.2) + stat_ellipse(type = "norm",color="#000000",data=f_gene_BRBclus3,size=0.2) + geom_point(aes(y=BRB, x=H3p3),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + xlab("H3.3") + ggtitle(pppplottitle)+ xlim(-1.0, 1.0) + ylim(-2.0, 2.0)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_nolabel_H3p3clus3_cutoff__H3p3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_H3p3clus3_cutoff__H3p3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3p3_BRB"),aes(x=-1.0,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3p3_BRB"),aes(x=-1.0,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3p3_BRB"),aes(x=-1.0,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_withcorr_H3p3clus3_cutoff__H3p3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3K4me3)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + stat_ellipse(type = "norm",color=density_color_high,size=0.2) + stat_ellipse(type = "norm",color="#000000",data=f_gene_BRBclus3,size=0.2) + geom_point(aes(y=BRB, x=H3K4me3),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + ggtitle(pppplottitle)+ xlim(-1.5, 1.5) + ylim(-2.0, 2.0)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_nolabel_H3p3clus3_cutoff__H3K4me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_H3p3clus3_cutoff__H3K4me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3K4me3_BRB"),aes(x=-1.5,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3K4me3_BRB"),aes(x=-1.5,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3K4me3_BRB"),aes(x=-1.5,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_withcorr_H3p3clus3_cutoff__H3K4me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3K27ac)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + stat_ellipse(type = "norm",color=density_color_high,size=0.2) + stat_ellipse(type = "norm",color="#000000",data=f_gene_BRBclus3,size=0.2) + geom_point(aes(y=BRB, x=H3K27ac),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + ggtitle(pppplottitle)+ xlim(-1.0, 1.0) + ylim(-2.0, 2.0)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_nolabel_H3p3clus3_cutoff__H3K27acvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_H3p3clus3_cutoff__H3K27acvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3K27ac_BRB"),aes(x=-1.0,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3K27ac_BRB"),aes(x=-1.0,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3K27ac_BRB"),aes(x=-1.0,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_withcorr_H3p3clus3_cutoff__H3K27acvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3K27me3)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + stat_ellipse(type = "norm",color=density_color_high,size=0.2) + stat_ellipse(type = "norm",color="#000000",data=f_gene_BRBclus3,size=0.2) + geom_point(aes(y=BRB, x=H3K27me3),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + ggtitle(pppplottitle)+ xlim(-1.0, 1.0) + ylim(-2.0, 2.0)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_nolabel_H3p3clus3_cutoff__H3K27me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_H3p3clus3_cutoff__H3K27me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3K27me3_BRB"),aes(x=-1.0,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3K27me3_BRB"),aes(x=-1.0,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3K27me3_BRB"),aes(x=-1.0,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_withcorr_H3p3clus3_cutoff__H3K27me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=ATAC)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + stat_ellipse(type = "norm",color=density_color_high,size=0.2) + stat_ellipse(type = "norm",color="#000000",data=f_gene_BRBclus3,size=0.2) + geom_point(aes(y=BRB, x=ATAC),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + ggtitle(pppplottitle)+ xlim(-0.4, 0.4) + ylim(-2.0, 2.0)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_nolabel_H3p3clus3_cutoff__ATACvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_H3p3clus3_cutoff__ATACvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="ATAC_BRB"),aes(x=-0.4,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="ATAC_BRB"),aes(x=-0.4,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="ATAC_BRB"),aes(x=-0.4,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/ellipse/log2FC_ellipse_withcorr_H3p3clus3_cutoff__ATACvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
For Check
#density_color_low <- "#ECE038"
density_color_low <- "#FFFFFF"
#density_color_high <- "#377EB8"
density_color_high <- "blue"
#density_color_low <- #FFFFFF"
#density_color_mid <- "yellow"
#density_color_high <- "red"
binsize <- 7
pppplottitle <- paste("log2 FC (Dox + vs -)\nBRB normalized count (Time, avg) > ",Set_cutoff,"\n H3.3 clus3: ",nrow(z_H3p3clus3)," genes\n Plot: ",H3p3clus3cutoff," genes\n BRB clus3: ",H3p3clus3cutoff_brbclus3," genes",sep="")
###
fcplot <-plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3p3)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density2d(aes(color=BRBDEGcluster),data=f_gene_BRBclus3,size=0.1,alpha = 0.5, bins=binsize) + geom_point(aes(y=BRB, x=H3p3, color=BRBDEGcluster,shape=shape),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + scale_shape_manual(values=c(21, 19)) + xlab("H3.3") + ggtitle(pppplottitle) + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_nolimit_H3p3clus3_cutoff__H3p3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + xlim(-1.0, 1.0) + ylim(-2.0, 2.0)
#fcplot
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_H3p3clus3_cutoff__H3p3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3p3_BRB"),aes(x=-1.0,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3p3_BRB"),aes(x=-1.0,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3p3_BRB"),aes(x=-1.0,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_withcorr_H3p3clus3_cutoff__H3p3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3K4me3)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density2d(aes(color=BRBDEGcluster),data=f_gene_BRBclus3,size=0.1,alpha = 0.5, bins=binsize) + geom_point(aes(y=BRB, x=H3K4me3, color=BRBDEGcluster,shape=shape),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + scale_shape_manual(values=c(21, 19)) + ggtitle(pppplottitle) + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_nolimit_H3p3clus3_cutoff__H3K4me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + xlim(-1.5, 1.5) + ylim(-2.0, 2.0)
#fcplot
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_H3p3clus3_cutoff__H3K4me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3K4me3_BRB"),aes(x=-1.5,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3K4me3_BRB"),aes(x=-1.5,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3K4me3_BRB"),aes(x=-1.5,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_withcorr_H3p3clus3_cutoff__H3K4me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3K27ac)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density2d(aes(color=BRBDEGcluster),data=f_gene_BRBclus3,size=0.1,alpha = 0.5, bins=binsize) + geom_point(aes(y=BRB, x=H3K27ac, color=BRBDEGcluster,shape=shape),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + scale_shape_manual(values=c(21, 19)) + ggtitle(pppplottitle) + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_nolimit_H3p3clus3_cutoff__H3K27acvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + xlim(-1.0, 1.0) + ylim(-2.0, 2.0)
#fcplot
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_H3p3clus3_cutoff__H3K27acvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3K27ac_BRB"),aes(x=-1.0,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3K27ac_BRB"),aes(x=-1.0,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3K27ac_BRB"),aes(x=-1.0,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_withcorr_H3p3clus3_cutoff__H3K27acvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=H3K27me3)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density2d(aes(color=BRBDEGcluster),data=f_gene_BRBclus3,size=0.1,alpha = 0.5, bins=binsize) + geom_point(aes(y=BRB, x=H3K27me3, color=BRBDEGcluster,shape=shape),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + scale_shape_manual(values=c(21, 19)) + ggtitle(pppplottitle) + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_nolimit_H3p3clus3_cutoff__H3K27me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + xlim(-1.0, 1.0) + ylim(-2.0, 2.0)
#fcplot
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_H3p3clus3_cutoff__H3K27me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="H3K27me3_BRB"),aes(x=-1.0,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="H3K27me3_BRB"),aes(x=-1.0,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="H3K27me3_BRB"),aes(x=-1.0,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_withcorr_H3p3clus3_cutoff__H3K27me3vsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(y=BRB, x=ATAC)) + facet_wrap(~time,ncol=1) + stat_density2d(aes(fill=..density..), geom = "raster",contour = FALSE) + scale_fill_gradient(low = density_color_low, high = density_color_high) + geom_abline(intercept=0,slope=0,colour="#000000",size=0.2) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density2d(aes(color=BRBDEGcluster),data=f_gene_BRBclus3,size=0.1,alpha = 0.5, bins=binsize) + geom_point(aes(y=BRB, x=ATAC, color=BRBDEGcluster,shape=shape),alpha = 0.6, size=1.0, data=f_gene_BRBclus3) + scale_color_manual(values = c("#000000"))+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=10),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + scale_shape_manual(values=c(21, 19)) + ggtitle(pppplottitle) + geom_text_repel(aes(label = label_text), segment.color = "#000000",segment.size = 0.1)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_nolimit_H3p3clus3_cutoff__ATACvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + xlim(-0.4, 0.4) + ylim(-2.0, 2.0)
#fcplot
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_H3p3clus3_cutoff__ATACvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot <- fcplot + geom_text_repel(data=filter(Cortest_H3p3clus3BRBclus3,Compare=="ATAC_BRB"),aes(x=-0.4,y=1.8,label=text), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(Cortest_H3p3clus3All,Compare=="ATAC_BRB"),aes(x=-0.4,y=2.0,label=text), color = density_color_high, segment.color = density_color_high,segment.size = 0.1,size = 1.8) + geom_text_repel(data=filter(summa__vsBRB_H3p3clus3_ALL_BRBclus3,Compare=="ATAC_BRB"),aes(x=-0.4,y=1.6,label=show), color = "#000000", segment.color = "#000000",segment.size = 0.1,size = 1.8)
ggsave(plot=fcplot,file="./log2FC/For_Check/log2FC_Check_withcorr_H3p3clus3_cutoff__ATACvsBRB.pdf", width = 3.5, height = 11, dpi = 360, limitsize = FALSE)
fcplot
20200817追加
#density_color_low <- "#ECE038"
density_color_low <- "#FFFFFF"
#density_color_high <- "#377EB8"
density_color_high <- "blue"
#density_color_low <- #FFFFFF"
#density_color_mid <- "yellow"
#density_color_high <- "red"
binsize <- 7
pppplottitle <- paste("log2 FC (Dox + vs -)\nBRB normalized count (Time, avg) > ",Set_cutoff,"\n H3.3 clus3: ",nrow(z_H3p3clus3)," genes\n Plot: ",H3p3clus3cutoff," genes\n BRB clus3: ",H3p3clus3cutoff_brbclus3," genes",sep="")
###
fcplot <-plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(x=BRB)) + facet_wrap(~time,ncol=1) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density(fill=density_color_high,color=density_color_high,alpha=0.5) + geom_density(fill="#000000",color="#000000",data=f_gene_BRBclus3,alpha=0.5) +theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=8),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank())+ xlim(-2.0, 2.0)
ggsave(plot=fcplot,file="./log2FC/histlog2FC__H3p3clus3_cutoff__BRB.pdf", width = 3.5, height = 3.5, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <-plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(x=H3p3)) + facet_wrap(~time,ncol=1) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density(fill=density_color_high,color=density_color_high,alpha=0.5) + geom_density(fill="#000000",color="#000000",data=f_gene_BRBclus3,alpha=0.5) +theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=8),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + xlab("H3.3")+ xlim(-1.0, 1.0)
ggsave(plot=fcplot,file="./log2FC/histlog2FC__H3p3clus3_cutoff__H3p3.pdf", width = 3.5, height = 3.5, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(x=H3K4me3)) + facet_wrap(~time,ncol=1) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density(fill=density_color_high,color=density_color_high,alpha=0.5) + geom_density(fill="#000000",color="#000000",data=f_gene_BRBclus3,alpha=0.5) +theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=8),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank())+ xlim(-1.5, 1.5)
ggsave(plot=fcplot,file="./log2FC/histlog2FC__H3p3clus3_cutoff__H3K4me3.pdf", width = 3.5, height = 3.5, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(x=H3K27ac)) + facet_wrap(~time,ncol=1) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density(fill=density_color_high,color=density_color_high,alpha=0.5) + geom_density(fill="#000000",color="#000000",data=f_gene_BRBclus3,alpha=0.5) +theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=8),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + xlim(-1.0, 1.0)
ggsave(plot=fcplot,file="./log2FC/histlog2FC__H3p3clus3_cutoff__H3K27ac.pdf", width = 3.5, height = 3.5, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(x=H3K27me3)) + facet_wrap(~time,ncol=1) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density(fill=density_color_high,color=density_color_high,alpha=0.5) + geom_density(fill="#000000",color="#000000",data=f_gene_BRBclus3,alpha=0.5) +theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=8),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + xlim(-1.0, 1.0)
ggsave(plot=fcplot,file="./log2FC/histlog2FC__H3p3clus3_cutoff__H3K27me3.pdf", width = 3.5, height = 3.5, dpi = 360, limitsize = FALSE)
fcplot
###
fcplot <- plot_all_FC_cutoff_H3p3clus3 %>% ggplot(aes(x=ATAC)) + facet_wrap(~time,ncol=1) + geom_vline(xintercept = 0,colour="#000000",size=0.2) + geom_density(fill=density_color_high,color=density_color_high,alpha=0.5) + geom_density(fill="#000000",color="#000000",data=f_gene_BRBclus3,alpha=0.5)+theme_bw() + theme(axis.title = element_text(size=15),axis.text = element_text(size=8),axis.text.x = element_text(hjust = 0.5,vjust=1.0), legend.position = "right", strip.text=element_text(size=15),strip.background = element_blank(),title = element_text(size=8),panel.grid=element_blank()) + xlim(-0.4, 0.4)
ggsave(plot=fcplot,file="./log2FC/histlog2FC__H3p3clus3_cutoff__ATAC.pdf", width = 3.5, height = 3.5, dpi = 360, limitsize = FALSE)
fcplot
2020.4.21, 7.17修正 ver
#20200421修正 ver
#20191206修正 ver
#z_heat_label_order_cluster6 <- z_heat_label_order_cluster %>% dplyr::select(ext_gene,heatmap_order,No,cluster_6) %>% mutate(heatmap_order=as.integer(heatmap_order),No=as.integer(No),cluster_6=as.integer(cluster_6))%>% arrange(heatmap_order) %>% left_join( dplyr::select(z_timedeg_s,ens_gene,ext_gene,biotype,chr))
#_____________#
## z_heat_label_order_cluster にクラスター番号が入っている
table_degcluster <- rrres_allH3p3 %>% filter(!is.na(cluster)) %>% arrange(cluster, ens_gene) %>% unique() %>% filter(!is.na(ens_gene))
degclusgene <- table_degcluster %>% group_by(cluster) %>% summarise(size=n()) %>% mutate(cluster=row_number())
table_degcluster <- table_degcluster %>% left_join(degclusgene %>% dplyr::select(cluster)) %>% arrange(cluster,ens_gene)
degclusgene
##### FDR setting ######
gofdr <- 0.1
#cluster_num <- 6
cluster_num <- nrow(degclusgene)
# 20191206修正
library(clusterProfiler)
library(org.Mm.eg.db)
folder_path <- "./H3p3allcluster/clusterProfile/"
#-------------#
file_path <- paste(folder_path, "GO_newcluster_BPfdr0p1_generatio",sep="")
filename_csv <- file_path
file_path <- paste(folder_path, "GO_newcluster_BPfdr0p1_generatio_cluster",sep="")
filename_list <- file_path
print(filename_list)
print(filename_csv)
#例 filename_list <- "./LRT/clusterProfile/H3mm18KO_mouseCTX_BRB0438_day5_2gunfdr0p2_kmeans_BPfdr0p1_generatio_cluster"
#例 filename_csv <- "./LRT/clusterProfile/H3mm18KO_mouseCTX_BRB0438_day5_2gunfdr0p2_kemans_BPfdr0p1_generatio"
#-------------#
cluster_list <- as.list(NA) #初期化
for (i in 1:cluster_num) {
pre_list <- as.list(NA)
pre_list <- table_degcluster %>% filter(cluster==as.integer(i)) %>% dplyr::select(ens_gene) %>% as.list()
names(pre_list) <- paste("ENSEMBL",as.character(i),sep="_")
if (i == 1) {
cluster_list <- pre_list
}
else cluster_list <- c(cluster_list, pre_list)
}
for (i in 1:cluster_num) {
print(paste(i, cluster_list[[i]] %>% tibble::enframe(name = NULL) %>% nrow(), sep=", "))
pre_ego_BP <- enrichGO(gene = cluster_list[[i]],
OrgDb = "org.Mm.eg.db",
keyType = 'ENSEMBL',
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = gofdr, qvalueCutoff = 1.0)
#20191211修正 pvalueCutoff = fdr
## pvalue < qvalue < p.adjust ##
# qvalueCutoff = 0.3 qvalueCutoff = 0.2 , qvalueCutoff = 1.0
#if (i == 1) {
# table_ego_BP <- data.frame(pre_ego_BP) %>% mutate(cluster=as.integer(i))
# # リスト型からデータフレームへ変換
#}
#else table_ego_BP <- table_ego_BP %>% bind_rows(data.frame(pre_ego_BP) %>% mutate(cluster=as.integer(i)))
if (i == 1) {
table_ego_BP <- data.frame(pre_ego_BP) %>% mutate(cluster=paste("cluster",as.character(i),sep="")) # リスト型からデータフレームへ変換
}
else table_ego_BP <- table_ego_BP %>% bind_rows(data.frame(pre_ego_BP) %>% mutate(cluster=paste("cluster",as.character(i),sep="")))
#---- plot ---#
BPplot <- dotplot(pre_ego_BP, showCategory=30, orderBy = "Count") #clusterProfile の機能で図を描く(191106修正) wrong orderBy parameter; set to default `orderBy = "x"`
print(BPplot)
ggsave(BPplot,file=paste(filename_list,as.character(i),".png",sep=""), width = 12, height = 12, dpi = 120)
ggsave(BPplot,file=paste(filename_list,as.character(i),".pdf",sep=""), width = 12, height = 12, dpi = 120)
}
print(table_ego_BP %>% group_by(cluster) %>% summarize())
#------#
# データはtable_ego_BPに。
#------------------------------------------------------#
# テーブルを保存
# table_ego_BP_3t3_LRT2 <- table_ego_BP
#
table_ego_BP1 <- table_ego_BP %>% mutate(cluster=factor(cluster,c("cluster1","cluster2","cluster3","cluster4","cluster5","cluster6"))) %>% arrange(cluster,desc(Count)) #191106(200415)
#table_ego_BP1 <- table_ego_BP %>% arrange(cluster,desc(Count)) %>% left_join(dplyr::select(degclusgene, cluster)) #191106(200415)
readr::write_csv(table_ego_BP1,paste(filename_csv,".csv",sep=""))
print(table_ego_BP %>% group_by(cluster) %>% summarize(cluster_3t3Dox_num = dplyr::n()))
# 先のテーブルのgeneIDをgene nameに置換する。(20191025)
tablego <- table_ego_BP1 %>% mutate(gene_name=geneID) %>% dplyr::select(-(qvalue))
for (i in 1:nrow(table_degcluster)) {
tablego <- tablego %>% mutate(gene_name=gsub(gene_name, pattern=table_degcluster$ens_gene[i], replacement=table_degcluster$ext_gene[i], ignore.case = TRUE))
}
#print(tablego)
#readr::write_csv(tablego,paste(filename_csv,"_genename.csv",sep=""))
#------------------------------------------------------#
readr::write_csv(tablego,paste(filename_csv,"_genename.csv",sep=""))
library(chromVAR)
library(motifmatchr)
library(SummarizedExperiment)
library(Matrix)
library(ggplot2)
library(BiocParallel)
library(BSgenome.Mmusculus.UCSC.mm10)
list_plotallFCcutoff_H3p3clus3 <- dplyr::select(plot_all_FC_cutoff_H3p3clus3, ens_gene, ext_gene, biotype, chr, H3p3cluster, BRBDEGcluster) %>% unique()
TSSregion_H3p3clus3 <- matome0_s %>% ungroup %>% dplyr::select(TSSstart,TSSend,ens_gene,score,strand,TSS,Start,End,position) %>% filter(ens_gene %in% list_plotallFCcutoff_H3p3clus3$ens_gene) %>% ungroup() %>% left_join(list_plotallFCcutoff_H3p3clus3)
Joining, by = "ens_gene"
#TSSPM10kb_H3p3clus3 <- TSSregion_H3p3clus3 %>% mutate(TSS_M10kb=TSS-10000,TSS_P10kb=TSS+10000) %>% dplyr::select(chr,TSS_M10kb,TSS_P10kb,ens_gene,score,strand,TSS, Start, End, position, ext_gene, biotype, H3p3cluster, BRBDEGcluster)
TSSPM10kb_H3p3clus3 <- TSSregion_H3p3clus3 %>% dplyr::select(chr,TSS,ens_gene, ext_gene,biotype, score, strand, H3p3cluster, BRBDEGcluster) %>% mutate(start=TSS-10000,end=TSS+10000) %>% dplyr::select(chr,start,end,ens_gene, score, strand,TSS, ext_gene,biotype, H3p3cluster, BRBDEGcluster) %>% filter(BRBDEGcluster=="3")
TSSPM10kb_H3p3clus3_GR <- makeGRangesFromDataFrame(TSSPM10kb_H3p3clus3, keep.extra.columns = TRUE)
seqlevelsStyle(TSSPM10kb_H3p3clus3_GR) <- "UCSC"
#####
peakfile <- "./Motif/TSSPM10kb_H3p3clus3.csv"
print(peakfile)
[1] "./Motif/TSSPM10kb_H3p3clus3.csv"
readr::write_csv(TSSPM10kb_H3p3clus3,peakfile)
head(TSSPM10kb_H3p3clus3)
nrow(TSSPM10kb_H3p3clus3)
[1] 55
def_bam_path <- "/home/guestA/o70578a/akuwakado/kuwakado/ChILSeq2/Komatsu_3T3_EGFP_H3mm18_Dox_chIl_0111NOVAseq/ChromVAR/ChromVAR_ChIL/H3K27ac_H3K27acpeak/deftable_ChromVAR_ChIL01100111_20200501_3T3_EGFP18_UI_DoxMinus_H3p3K27acK4Kme327me3.txt"
def_H3K27ac_bam <- readr::read_tsv(file =def_bam_path) %>% filter(seq=="H3K27ac")
Parsed with column specification:
cols(
file = col_character(),
multicov_No = col_double(),
sample = col_character(),
group = col_character(),
time = col_character(),
type = col_character(),
seq = col_character(),
rep = col_double()
)
bamfiles <- def_H3K27ac_bam$file #bamfiles <- def_bam$peakcall_bam
fragment_counts <- getCounts(def_H3K27ac_bam$file, TSSPM10kb_H3p3clus3_GR,
paired = FALSE, # ChILはペアではない。
by_rg = FALSE,
colData = DataFrame(Cell_Type = def_H3K27ac_bam$group,sample_name = def_H3K27ac_bam$sample))
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_UI_DoxMinus_H3K27ac_Rep01.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_UI_DoxMinus_H3K27ac_Rep02.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_UI_DoxPlus_H3K27ac_Rep01.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_UI_DoxPlus_H3K27ac_Rep02.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep01.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep02.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_0h_DoxPlus_H3K27ac_Rep01.bam
Reading in file: /home/guestA/n70275a/kTanaka/0110NOVAseq/mapped/20200501_3T3_EGFP18_0h_DoxPlus_H3K27ac_Rep02.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_24h_DoxMinus_H3K27ac_Rep01_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_24h_DoxMinus_H3K27ac_Rep02_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_24h_DoxPlus_H3K27ac_Rep01_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_24h_DoxPlus_H3K27ac_Rep02_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep01_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep02_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep01_runMerged.bam
Reading in file: /home/guestA/n70275a/kTanaka/0111NOVAseq/mapped/3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep02_runMerged.bam
length(fragment_counts)
[1] 55
print("---- fragment_counts ----")
[1] "---- fragment_counts ----"
slot(fragment_counts, "rowRanges") #fragment_counts @rowRanges
GRanges object with 55 ranges and 7 metadata columns:
seqnames ranges strand | ens_gene score TSS ext_gene biotype H3p3cluster
<Rle> <IRanges> <Rle> | <character> <character> <integer> <character> <character> <integer>
[1] chr1 51279126-51299126 + | ENSMUSG00000045954 . 51289126 Cavin2 protein_coding 3
[2] chr1 75350329-75370329 + | ENSMUSG00000026208 . 75360329 Des protein_coding 3
[3] chr1 134279989-134299989 + | ENSMUSG00000026459 . 134289989 Myog protein_coding 3
[4] chr1 134322928-134342928 - | ENSMUSG00000026458 . 134332928 Ppfia4 protein_coding 3
[5] chr1 135365237-135385237 - | ENSMUSG00000041889 . 135375237 Shisa4 protein_coding 3
... ... ... ... . ... ... ... ... ... ...
[51] chr19 6986117-7006117 - | ENSMUSG00000037349 . 6996117 Nudt22 protein_coding 3
[52] chr19 34245590-34265590 - | ENSMUSG00000035783 . 34255590 Acta2 protein_coding 3
[53] chr19 40503779-40523779 - | ENSMUSG00000025006 . 40513779 Sorbs1 protein_coding 3
[54] chr19 42026000-42046000 + | ENSMUSG00000025172 . 42036000 Ankrd2 protein_coding 3
[55] chrX 167199315-167219315 - | ENSMUSG00000049775 . 167209315 Tmsb4x protein_coding 3
BRBDEGcluster
<factor>
[1] 3
[2] 3
[3] 3
[4] 3
[5] 3
... ...
[51] 3
[52] 3
[53] 3
[54] 3
[55] 3
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
slot(fragment_counts, "colData")
DataFrame with 16 rows and 3 columns
Cell_Type sample_name depth
<character> <character> <numeric>
20200501_3T3_EGFP18_UI_DoxMinus_H3K27ac_Rep01.bam H3K27ac_UI_DoxMinus H3K27ac_UI_DoxMinus_1 2655505
20200501_3T3_EGFP18_UI_DoxMinus_H3K27ac_Rep02.bam H3K27ac_UI_DoxMinus H3K27ac_UI_DoxMinus_2 3467576
20200501_3T3_EGFP18_UI_DoxPlus_H3K27ac_Rep01.bam H3K27ac_UI_DoxPlus H3K27ac_UI_DoxPlus_1 3198781
20200501_3T3_EGFP18_UI_DoxPlus_H3K27ac_Rep02.bam H3K27ac_UI_DoxPlus H3K27ac_UI_DoxPlus_2 3038383
20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep01.bam H3K27ac_0h_DoxMinus H3K27ac_0h_DoxMinus_1 1959727
... ... ... ...
3T3_EGFP18_24h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_24h_DoxPlus H3K27ac_24h_DoxPlus_2 1527310
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_1 1101851
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_2 1054200
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_1 1241266
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_2 889154
slot(fragment_counts, "NAMES")
NULL
slot(fragment_counts, "metadata")
list()
slot(fragment_counts, "elementMetadata")
DataFrame with 55 rows and 0 columns
slot(fragment_counts, "assays")
An object of class "SimpleAssays"
Slot "data":
List of length 1
names(1): counts
print("--------- save ----------------")
[1] "--------- save ----------------"
#-- save fragment count ----#
#--------#
type_depth <- slot(fragment_counts, "colData") %>% as.data.frame() # 範囲
fffile <- sub(".csv","_typedepth.csv",peakfile)
print(fffile)
[1] "./Motif/TSSPM10kb_H3p3clus3_typedepth.csv"
type_depth %>% readr::write_csv(fffile)
#--------#
f_c_range <- fragment_counts @rowRanges %>% as.data.frame() # 範囲
fc_range <- f_c_range %>% mutate(seqnames1=seqnames,start1=start,end1=end) %>% unite(sten,c(start1,end1),sep="-") %>% unite(range,c(seqnames1,sten),sep=":")
f_c_count <- fragment_counts @assays @data$counts %>% as.matrix() %>% as.data.frame() # カウント
#---#
f_c_range_count <- cbind(fc_range, f_c_count)
nrow(f_c_range_count)
[1] 55
fffile <- sub(".csv","_fragcounts.csv",peakfile)
print(fffile)
[1] "./Motif/TSSPM10kb_H3p3clus3_fragcounts.csv"
f_c_range_count %>% readr::write_csv(fffile)
#---------------------------#
register(SerialParam())
fragment_counts_bias <- addGCBias(fragment_counts, genome = BSgenome.Mmusculus.UCSC.mm10) #ここでこけないように、Chrは確定されているものに設定。
#+++++++++++++++++++++++++++++++++++++++++++++++#
length(fragment_counts_bias)
[1] 55
print("---- fragment_counts_bias ----")
[1] "---- fragment_counts_bias ----"
print("== rowRanges ==")
[1] "== rowRanges =="
slot(fragment_counts_bias, "rowRanges")
GRanges object with 55 ranges and 8 metadata columns:
seqnames ranges strand | ens_gene score TSS ext_gene biotype H3p3cluster
<Rle> <IRanges> <Rle> | <character> <character> <integer> <character> <character> <integer>
[1] chr1 51279126-51299126 + | ENSMUSG00000045954 . 51289126 Cavin2 protein_coding 3
[2] chr1 75350329-75370329 + | ENSMUSG00000026208 . 75360329 Des protein_coding 3
[3] chr1 134279989-134299989 + | ENSMUSG00000026459 . 134289989 Myog protein_coding 3
[4] chr1 134322928-134342928 - | ENSMUSG00000026458 . 134332928 Ppfia4 protein_coding 3
[5] chr1 135365237-135385237 - | ENSMUSG00000041889 . 135375237 Shisa4 protein_coding 3
... ... ... ... . ... ... ... ... ... ...
[51] chr19 6986117-7006117 - | ENSMUSG00000037349 . 6996117 Nudt22 protein_coding 3
[52] chr19 34245590-34265590 - | ENSMUSG00000035783 . 34255590 Acta2 protein_coding 3
[53] chr19 40503779-40523779 - | ENSMUSG00000025006 . 40513779 Sorbs1 protein_coding 3
[54] chr19 42026000-42046000 + | ENSMUSG00000025172 . 42036000 Ankrd2 protein_coding 3
[55] chrX 167199315-167219315 - | ENSMUSG00000049775 . 167209315 Tmsb4x protein_coding 3
BRBDEGcluster bias
<factor> <numeric>
[1] 3 0.418129
[2] 3 0.492075
[3] 3 0.498425
[4] 3 0.521024
[5] 3 0.501925
... ... ...
[51] 3 0.535723
[52] 3 0.421279
[53] 3 0.459677
[54] 3 0.495675
[55] 3 0.429779
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
print("== colData ==")
[1] "== colData =="
slot(fragment_counts_bias, "colData")
DataFrame with 16 rows and 3 columns
Cell_Type sample_name depth
<character> <character> <numeric>
20200501_3T3_EGFP18_UI_DoxMinus_H3K27ac_Rep01.bam H3K27ac_UI_DoxMinus H3K27ac_UI_DoxMinus_1 2655505
20200501_3T3_EGFP18_UI_DoxMinus_H3K27ac_Rep02.bam H3K27ac_UI_DoxMinus H3K27ac_UI_DoxMinus_2 3467576
20200501_3T3_EGFP18_UI_DoxPlus_H3K27ac_Rep01.bam H3K27ac_UI_DoxPlus H3K27ac_UI_DoxPlus_1 3198781
20200501_3T3_EGFP18_UI_DoxPlus_H3K27ac_Rep02.bam H3K27ac_UI_DoxPlus H3K27ac_UI_DoxPlus_2 3038383
20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep01.bam H3K27ac_0h_DoxMinus H3K27ac_0h_DoxMinus_1 1959727
... ... ... ...
3T3_EGFP18_24h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_24h_DoxPlus H3K27ac_24h_DoxPlus_2 1527310
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_1 1101851
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_2 1054200
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_1 1241266
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_2 889154
print("== NAMES ==")
[1] "== NAMES =="
slot(fragment_counts_bias, "NAMES")
NULL
print("== metadata ==")
[1] "== metadata =="
slot(fragment_counts_bias, "metadata")
list()
print("== elementMetadata ==")
[1] "== elementMetadata =="
slot(fragment_counts_bias, "elementMetadata")
DataFrame with 55 rows and 0 columns
print("== assays ==")
[1] "== assays =="
slot(fragment_counts_bias, "assays")
An object of class "SimpleAssays"
Slot "data":
List of length 1
names(1): counts
print("--------- save ----------------")
[1] "--------- save ----------------"
#-- save fragment count bias----#
f_c_range_bias <- fragment_counts_bias @rowRanges %>% as.data.frame() # 範囲
fc_range_bias <- f_c_range_bias %>% mutate(seqnames1=seqnames,start1=start,end1=end) %>% unite(sten,c(start1,end1),sep="-") %>% unite(range,c(seqnames1,sten),sep=":")
f_c_count_bias <- fragment_counts_bias @assays @data$counts %>% as.matrix() %>% as.data.frame() # カウント
#---#
f_c_range_count_bias <- cbind(fc_range_bias, f_c_count_bias)
nrow(f_c_range_count_bias)
[1] 55
fffile <- sub(".csv","_fragcounts_bias.csv",peakfile)
print(fffile)
[1] "./Motif/TSSPM10kb_H3p3clus3_fragcounts_bias.csv"
f_c_range_count_bias %>% readr::write_csv(fffile)
#---------------------------#
#counts_filtered <- filterSamples(fragment_counts_bias, min_depth = 1500, min_in_peaks = 0.15, shiny = FALSE)
counts_filtered_pre <- filterSamples(fragment_counts_bias, shiny = FALSE)
min_in_peaks set to 0.002
min_depth set to 161032.75
#++++++++++++++++++++++++++++++++++++++++++++++#
# If unspecified, min_in_peaks and min_depth cutoffs will be estimated based on data. min_in_peaks is set to 0.5 times the median proportion of fragments in peaks. min_depth is set to the maximum of 500 or 10 median library size.
#
# min_in_peaks: minimum fraction of samples within peaks
# min_depth: minimum library size
# shiny: make shiny gadget?
# ix_return: return indices of sample to keep instead of subsetted counts object
#++++++++++++++++++++++++++++++++++++++++++++++#
length(counts_filtered_pre)
[1] 55
print("---- counts_filtered (pre) ----")
[1] "---- counts_filtered (pre) ----"
print("== rowRanges ==")
[1] "== rowRanges =="
slot(counts_filtered_pre, "rowRanges")
GRanges object with 55 ranges and 8 metadata columns:
seqnames ranges strand | ens_gene score TSS ext_gene biotype H3p3cluster
<Rle> <IRanges> <Rle> | <character> <character> <integer> <character> <character> <integer>
[1] chr1 51279126-51299126 + | ENSMUSG00000045954 . 51289126 Cavin2 protein_coding 3
[2] chr1 75350329-75370329 + | ENSMUSG00000026208 . 75360329 Des protein_coding 3
[3] chr1 134279989-134299989 + | ENSMUSG00000026459 . 134289989 Myog protein_coding 3
[4] chr1 134322928-134342928 - | ENSMUSG00000026458 . 134332928 Ppfia4 protein_coding 3
[5] chr1 135365237-135385237 - | ENSMUSG00000041889 . 135375237 Shisa4 protein_coding 3
... ... ... ... . ... ... ... ... ... ...
[51] chr19 6986117-7006117 - | ENSMUSG00000037349 . 6996117 Nudt22 protein_coding 3
[52] chr19 34245590-34265590 - | ENSMUSG00000035783 . 34255590 Acta2 protein_coding 3
[53] chr19 40503779-40523779 - | ENSMUSG00000025006 . 40513779 Sorbs1 protein_coding 3
[54] chr19 42026000-42046000 + | ENSMUSG00000025172 . 42036000 Ankrd2 protein_coding 3
[55] chrX 167199315-167219315 - | ENSMUSG00000049775 . 167209315 Tmsb4x protein_coding 3
BRBDEGcluster bias
<factor> <numeric>
[1] 3 0.418129
[2] 3 0.492075
[3] 3 0.498425
[4] 3 0.521024
[5] 3 0.501925
... ... ...
[51] 3 0.535723
[52] 3 0.421279
[53] 3 0.459677
[54] 3 0.495675
[55] 3 0.429779
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
print("== colData ==")
[1] "== colData =="
slot(counts_filtered_pre, "colData")
DataFrame with 12 rows and 3 columns
Cell_Type sample_name depth
<character> <character> <numeric>
20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep01.bam H3K27ac_0h_DoxMinus H3K27ac_0h_DoxMinus_1 1959727
20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep02.bam H3K27ac_0h_DoxMinus H3K27ac_0h_DoxMinus_2 1693345
20200501_3T3_EGFP18_0h_DoxPlus_H3K27ac_Rep01.bam H3K27ac_0h_DoxPlus H3K27ac_0h_DoxPlus_1 3424694
20200501_3T3_EGFP18_0h_DoxPlus_H3K27ac_Rep02.bam H3K27ac_0h_DoxPlus H3K27ac_0h_DoxPlus_2 2490017
3T3_EGFP18_24h_DoxMinus_H3K27ac_Rep01_runMerged.bam H3K27ac_24h_DoxMinus H3K27ac_24h_DoxMinus_1 1030885
... ... ... ...
3T3_EGFP18_24h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_24h_DoxPlus H3K27ac_24h_DoxPlus_2 1527310
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_1 1101851
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_2 1054200
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_1 1241266
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_2 889154
print("== NAMES ==")
[1] "== NAMES =="
slot(counts_filtered_pre, "NAMES")
NULL
print("== metadata ==")
[1] "== metadata =="
slot(counts_filtered_pre, "metadata")
list()
print("== elementMetadata ==")
[1] "== elementMetadata =="
slot(counts_filtered_pre, "elementMetadata")
DataFrame with 55 rows and 0 columns
print("== assays ==")
[1] "== assays =="
slot(counts_filtered_pre, "assays")
An object of class "SimpleAssays"
Slot "data":
List of length 1
names(1): counts
print("-------------------------")
[1] "-------------------------"
counts_filtered <- filterPeaks(counts_filtered_pre,non_overlapping=TRUE)
#++++++++++++++++++++++++++++++++++++++++++++++#
# if non_overlapping is set to true, when peaks overlap the overlapping peak with lower counts is removed
#
# min_fragments_per_peak: minimum number of fragmints in peaks across all samples
# non_overlapping: reduce peak set to non-overlapping peaks, see details
# ix_return: return indices of peaks to keep instead of subsetted counts object
#++++++++++++++++++++++++++++++++++++++++++++++#
length(counts_filtered)
[1] 55
print("---- counts_filtered ----")
[1] "---- counts_filtered ----"
slot(counts_filtered, "rowRanges")
GRanges object with 55 ranges and 8 metadata columns:
seqnames ranges strand | ens_gene score TSS ext_gene biotype H3p3cluster
<Rle> <IRanges> <Rle> | <character> <character> <integer> <character> <character> <integer>
[1] chr1 51279126-51299126 + | ENSMUSG00000045954 . 51289126 Cavin2 protein_coding 3
[2] chr1 75350329-75370329 + | ENSMUSG00000026208 . 75360329 Des protein_coding 3
[3] chr1 134279989-134299989 + | ENSMUSG00000026459 . 134289989 Myog protein_coding 3
[4] chr1 134322928-134342928 - | ENSMUSG00000026458 . 134332928 Ppfia4 protein_coding 3
[5] chr1 135365237-135385237 - | ENSMUSG00000041889 . 135375237 Shisa4 protein_coding 3
... ... ... ... . ... ... ... ... ... ...
[51] chr19 6986117-7006117 - | ENSMUSG00000037349 . 6996117 Nudt22 protein_coding 3
[52] chr19 34245590-34265590 - | ENSMUSG00000035783 . 34255590 Acta2 protein_coding 3
[53] chr19 40503779-40523779 - | ENSMUSG00000025006 . 40513779 Sorbs1 protein_coding 3
[54] chr19 42026000-42046000 + | ENSMUSG00000025172 . 42036000 Ankrd2 protein_coding 3
[55] chrX 167199315-167219315 - | ENSMUSG00000049775 . 167209315 Tmsb4x protein_coding 3
BRBDEGcluster bias
<factor> <numeric>
[1] 3 0.418129
[2] 3 0.492075
[3] 3 0.498425
[4] 3 0.521024
[5] 3 0.501925
... ... ...
[51] 3 0.535723
[52] 3 0.421279
[53] 3 0.459677
[54] 3 0.495675
[55] 3 0.429779
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
slot(counts_filtered, "colData")
DataFrame with 12 rows and 3 columns
Cell_Type sample_name depth
<character> <character> <numeric>
20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep01.bam H3K27ac_0h_DoxMinus H3K27ac_0h_DoxMinus_1 1959727
20200501_3T3_EGFP18_0h_DoxMinus_H3K27ac_Rep02.bam H3K27ac_0h_DoxMinus H3K27ac_0h_DoxMinus_2 1693345
20200501_3T3_EGFP18_0h_DoxPlus_H3K27ac_Rep01.bam H3K27ac_0h_DoxPlus H3K27ac_0h_DoxPlus_1 3424694
20200501_3T3_EGFP18_0h_DoxPlus_H3K27ac_Rep02.bam H3K27ac_0h_DoxPlus H3K27ac_0h_DoxPlus_2 2490017
3T3_EGFP18_24h_DoxMinus_H3K27ac_Rep01_runMerged.bam H3K27ac_24h_DoxMinus H3K27ac_24h_DoxMinus_1 1030885
... ... ... ...
3T3_EGFP18_24h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_24h_DoxPlus H3K27ac_24h_DoxPlus_2 1527310
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_1 1101851
3T3_EGFP18_48h_DoxMinus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxMinus H3K27ac_48h_DoxMinus_2 1054200
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep01_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_1 1241266
3T3_EGFP18_48h_DoxPlus_H3K27ac_Rep02_runMerged.bam H3K27ac_48h_DoxPlus H3K27ac_48h_DoxPlus_2 889154
slot(counts_filtered, "NAMES")
NULL
slot(counts_filtered, "metadata")
list()
slot(counts_filtered, "elementMetadata")
DataFrame with 55 rows and 0 columns
slot(counts_filtered, "assays")
An object of class "SimpleAssays"
Slot "data":
List of length 1
names(1): counts
print("-------------------------")
[1] "-------------------------"
#-- save counts_filtered ----#
f_c_range_countsfil <- counts_filtered @rowRanges %>% as.data.frame() # 範囲
fc_range_countsfil <- f_c_range_countsfil %>% mutate(seqnames1=seqnames,start1=start,end1=end) %>% unite(sten,c(start1,end1),sep="-") %>% unite(range,c(seqnames1,sten),sep=":")
f_c_count_countsfil <- counts_filtered @assays @data$counts %>% as.matrix() %>% as.data.frame() # カウント
#---#
f_c_range_count_countsfil <- cbind(fc_range_countsfil, f_c_count_countsfil)
nrow(f_c_range_count_countsfil)
[1] 55
fffile <- sub(".csv","_countsfilter.csv",peakfile)
print(fffile)
[1] "./Motif/TSSPM10kb_H3p3clus3_countsfilter.csv"
f_c_range_count_countsfil %>% readr::write_csv(fffile)
#---------------------------#
Raw deviations for background peaks & Bias corrected deviations and Z-scores
length(counts_filtered)
[1] 55
motifs <- chromVAR::getJasparMotifs(species = "Mus musculus", collection = "CORE") #OK
motif_ix <- matchMotifs(motifs, counts_filtered, genome = BSgenome.Mmusculus.UCSC.mm10)
#motifMatches(motif_ix) # Extract matches matrix from SummarizedExperiment result
dev <- computeDeviations(object = counts_filtered, annotations = motif_ix)
.doLoadActions(where, attach) でエラー:
error in load action .__A__.1 for package nabor: loadModule(module, NULL, env = where, loadNow = TRUE): Unable to load module "class_WKNNF": サイズ 1040604.2 Gb のベクトルを割り当てることができません
ChromVARで エラーのためここで終了(20200819)
(20200617追加) https://greenleaflab.github.io/chromVAR/reference/differentialDeviations.html
data(mini_dev, package = “chromVAR”) difdev <- differentialDeviations(mini_dev, “Cell_Type”) differentialDeviations(object, groups, alternative = c(“two.sided”, “less”,“greater”), parametric = TRUE) dev@colData
difdev <- differentialDeviations(dev, "Cell_Type", alternative = c("two.sided"))
is(object, "chromVARDeviations") でエラー:
オブジェクト 'dev' がありません
%>% filter(!motif_name %in% c(“Myod1”,“Myog”,“Tcf12”,“Tcf21”,“Ascl2”)) %>% filter(!motif_name %in% c(“FOS::JUN”,“Nfe2l2”,“Bach1::Mafk”)) %>% filter(!motif_name %in% c(“RUNX1”,“Myb”)) %>% filter(!motif_name %in% c(“Bcl6”,“Klf12”,“Klf4”,“Klf1”,“Gata4”,“Gata1”,“Rfx1”,“Spz1”,“Myc”,“Atoh1”))
%>% filter(motif_name %in% c(“Myod1”,“Myog”,“Tcf12”,“Tcf21”,“Ascl2”)) %>% filter(motif_name %in% c(“FOS::JUN”,“Nfe2l2”,“Bach1::Mafk”)) %>% filter(motif_name %in% c(“RUNX1”,“Myb”)) %>% filter(motif_name %in% c(“Bcl6”,“Klf12”,“Klf4”,“Klf1”,“Gata4”,“Gata1”,“Rfx1”,“Spz1”,“Myc”,“Atoh1”))
#+++++++++++++++++++++++++++++++++++++++++++++++#
print("---- motif_ix ----------------")
[1] "---- motif_ix ----------------"
motif_id <- motif_ix @colData %>% as_tibble(rownames = "motif_ID") %>% dplyr::rename(motif_name=name) #rename modif 20200616
print("-- motifMatches(motif_ix) --")
[1] "-- motifMatches(motif_ix) --"
## motifMatches が . or | で 入っている
motifM_table <- motifMatches(motif_ix) %>% as.matrix() %>% as.data.frame() #motifM_table <- motifMatches(motif_ix) %>% as.matrix() %>% as.data.frame() as_tibble()
#motifM_table %>% dplyr::select(1:3)
#colnames(motifM_table) #colnames(motifMatches(motif_ix))
ncol(motifM_table) #ncol(motifMatches(motif_ix))
[1] 128
motifM_table_range <- cbind(fc_range_countsfil, motifM_table)
ncol(motifM_table_range)
[1] 142
#motifM_table_range <- cbind(fc_range_bias %>% dplyr::select(range) ,motifM_table)
print("---- save motif ----")
[1] "---- save motif ----"
#-- save motif id ----#
fffile <- paste("./bed_allpeak/",sub(".bed","_maxs_allpeak_deg_500_motifID.csv",basename(peakfile)),sep="")
print(fffile)
[1] "./bed_allpeak/TSSPM10kb_H3p3clus3.csv"
motif_id %>% readr::write_csv(fffile)
ファイル './bed_allpeak/TSSPM10kb_H3p3clus3.csv' を開くことができません: そのようなファイルやディレクトリはありません open.connection(path, "wb") でエラー:
コネクションを開くことができません